Symbiont ecology and its influence on the physiology and isotopic composition of a reef coral across space and time
CB Wall cbwall@hawaii.edu
9/19/2018
require("knitr")
knitr::opts_chunk$set(warning=FALSE, message=FALSE, fig.align='center')
library(ggplot2)
library(broom)
library(tidyverse)
library(plyr)
library(dplyr)
library(plotrix)
library(lubridate)
library(effects)
library(RCurl)
library(chron)
library(lubridate)
library(lme4)
library(emmeans)
library(multcomp)
library(devtools)
library(logihist)
library(popbio)
library(car)
library(lmerTest)
library(cowplot)
library(sjPlot)
library(sjmisc)Project background
Environmental
Sea level correction
Corals were collected across in summer and winter 2016 across several different days with each period. To correct for differences in the depth of coral colonies we used NOAA tide data to correct all depth measurements to Mean Sea Level (MSL). This approach was done in Innis et al. 2018 on a larger collection of the coral samples–which have been subset in the current study.
#########################
# Sea level correction
#########################
rm(list=ls())
#### attach data files
data<-read.csv("data/mastersheet_PanKBAY.csv") # master file
qPCR.Innis<-read.csv("data/PanKBay_summer_qPCR.csv") # qPCR from summer (Innis et al. 2018)
JuneTide=read.csv("data/environmental/sea level/Station_1612480_tide_ht_20160601-20160630.csv")
JulyTide=read.csv("data/environmental/sea level/Station_1612480_tide_ht_20160701-20160731.csv")
AugustTide=read.csv("data/environmental/sea level/Station_1612480_tide_ht_20160801-20160831.csv")
DecemberTide=read.csv("data/environmental/sea level/Station_1612480_tide_ht_20161201-20161222.csv")
Tide<-rbind(JuneTide, JulyTide, AugustTide, DecemberTide) # all MSL in meters
Tide$Time <- as.POSIXct(Tide$TimeUTC, format="%Y-%m-%d %H:%M:%S", tz="UTC")
attributes(Tide$Time)$tzone <- "Pacific/Honolulu"
# use attr(as.POSIXlt(Tide$Time),"tzone") to confirm in PST
Tide<-Tide[, c(-5:-8)] #remove unnecessary columns
data$date.time<- as.POSIXct(paste(as.character(data$Date), as.character(data$Time.of.collection)),
format="%m/%d/%y %H:%M", tz="Pacific/Honolulu") # make sure time in HST for master
data$Time.r<-round_date(data$date.time,unit="6 minutes")
Tide$Time.r <- Tide$Time
data<-merge(data, Tide, by="Time.r", all.x=T)
data$newDepth <- data$Depth..m-data$TideHT # new depth in METERSLight parameters
DLI 2m all sites
This figure shows the daily integrated light (DLI mol photon m-2 d-1) at long-term light loggers at 2m depth near the locations where corals were collected in northern and southern sectors of Kāne’ohe Bay near shore (east) further off shore (west).
# DLI
##### all DLI for graphs
files <- list.files(path="data/environmental/temp and light/Jun_DecPAR/all DLI", pattern = "csv$", full.names = T)
tables <- lapply(files, read.csv, header = TRUE)
DLI<-do.call(rbind, tables)
DLI=DLI[-1] # remove junk column
DLI$Date<-as.Date(DLI$Date) # fix date
DLI<-DLI[order(DLI$Date),] # order by Date
# make a date sequence for entire study
all.date<-as.data.frame(seq(as.Date("2016-06-10"), as.Date("2017-01-12"), "days"))
colnames(all.date)="Date"
# merge the DLI data and the date sequence to make a complete df through time
DLI.3site<-merge(all.date, DLI, by="Date", all.x=T)
R10<-read.csv("data/environmental/temp and light/Jun_DecPAR/all DLI/Reef 10/DLI.Rf102016.csv")
R10<-R10[-1]; R10$Date<-as.Date(R10$Date)
DLI.4site<-merge(DLI.3site, R10, by="Date", all.x=T)
DLI.4site$month <- months(as.Date(DLI.4site$Date)) # makes a month column
DLI.4site<-DLI.4site[, c(1,6, 2:5)]
# write.csv(DLI.4site, "data/environmental/temp and light/Jun_DecPAR/All.DLI.csv")
##### Figure
All.DLI<-DLI.4site
reefcols=c("mediumseagreen", "dodgerblue", "salmon", "orchid")
par(mar=c(2,3.6,1,1.3), mgp=c(2,0.5,0))
layout(matrix(c(1,2,3,4), 2,2, byrow=TRUE))
plot(Rf44.mid.DLI~Date, type="l", All.DLI, ylab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), ylim=c(0, 30), xaxt="n", xlab="Time", col=reefcols[1])
legend("topright", lty=1, col=reefcols[1], legend="Northwest", lwd=1.5, bty="n", cex=0.8,
x.intersp = 0.3, inset=c(0.02, -0.02))
axis.Date(1, at=seq(min(All.DLI$Date), max(All.DLI$Date), by="2 mo"), format="%b '%y")
plot(Rf42.mid.DLI~Date, type="l", All.DLI, ylab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))), ylim=c(0, 30), xaxt="n", xlab="Time", col=reefcols[2])
legend("topright", lty=1, col=reefcols[2], legend="Northeast", lwd=1.5, bty="n", cex=0.8,
x.intersp = 0.3, inset=c(0.02, -0.02))
axis.Date(1, at=seq(min(All.DLI$Date), max(All.DLI$Date), by="2 mo"), format="%b '%y")
plot(Rf10.mid.DLI~Date, type="l", All.DLI, ylab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1, ")", sep=""))), ylim=c(0, 30), xaxt="n", xlab="Time", col=reefcols[3])
legend("topright", lty=1, col=reefcols[3], legend="Southwest", lwd=1.5, bty="n", cex=0.8,
x.intersp = 0.3, inset=c(0.02, -0.02))
axis.Date(1, at=seq(min(All.DLI$Date), max(All.DLI$Date), by="2 mo"), format="%b '%y")
plot(HIMB.mid.DLI~Date, type="l", All.DLI, ylab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1, ")", sep=""))), ylim=c(0, 30), xaxt="n", xlab="Time", col=reefcols[4])
legend("topright", lty=1, col=reefcols[4], legend="Southeast", lwd=1.5, bty="n", cex=0.8,
x.intersp = 0.3, inset=c(0.02, -0.02))
axis.Date(1, at=seq(min(All.DLI$Date), max(All.DLI$Date), by="2 mo"), format="%b '%y")
Figure DLI at the four reef areas from summer to winter 2016 at 2m depths. Data at 2m was recorded from loggers
dev.copy(pdf, "figures/environmental/All.DLI.pdf", width=6.5, height=4)
dev.off()Model results testing differences in DLI at 2 m depth at all sites.
mod.light<-DLI.4site
Rf44DLI.mod<-mod.light[1:3]; Rf44DLI.mod$Location<-"Northwest";
colnames(Rf44DLI.mod)<-c("Date", "month", "DLI", "Location")
Rf42DLI.mod<-mod.light[,c(1:2,4)]; Rf42DLI.mod$Location<-"Northeast";
colnames(Rf42DLI.mod)<-c("Date", "month", "DLI", "Location")
HIMBDLI.mod<-mod.light[,c(1:2,5)]; HIMBDLI.mod$Location<-"Southeast";
colnames(HIMBDLI.mod)<-c("Date", "month", "DLI", "Location")
Rf10DLI.mod<-mod.light[,c(1:2,6)]; Rf10DLI.mod$Location<-"Southwest";
colnames(Rf10DLI.mod)<-c("Date", "month", "DLI", "Location")
All.DLI.mod<-rbind(Rf44DLI.mod, Rf42DLI.mod, HIMBDLI.mod, Rf10DLI.mod)
All.DLI.mod$Location<-as.factor(All.DLI.mod$Location)
All.DLI.mod$month<-as.factor(All.DLI.mod$month)
All.DLI.mod$Season<-ifelse(All.DLI.mod$month== "June" |
All.DLI.mod$month== "July" |
All.DLI.mod$month== "August" |
All.DLI.mod$month== "September", "dry season", "wet season")
All.DLI.mod$Season<-as.factor(All.DLI.mod$Season)
mod<-lmer(DLI~Location*Season + (1|Date), data=All.DLI.mod); anova(mod, type=2)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Location 4378.8 1459.58 3 529.66 134.674 < 2.2e-16 ***
## Season 1040.9 1040.91 1 210.48 96.043 < 2.2e-16 ***
## Location:Season 490.9 163.62 3 531.14 15.097 1.924e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(allEffects(mod))
posthoc<-emmeans(mod, ~Location|Season)
CLD(posthoc, Letters=letters)## Season = dry season:
## Location emmean SE df lower.CL upper.CL .group
## Southwest 8.149681 0.4014748 504.55 7.360912 8.938449 a
## Northwest 10.829854 0.4450350 590.91 9.955811 11.703896 b
## Southeast 14.748648 0.4016675 503.81 13.959498 15.537797 c
## Northeast 14.997825 0.4601858 615.62 14.094100 15.901549 c
##
## Season = wet season:
## Location emmean SE df lower.CL upper.CL .group
## Southwest 3.901772 0.4473755 539.33 3.022959 4.780584 a
## Southeast 8.607241 0.4662620 551.22 7.691373 9.523109 b
## Northwest 9.128302 0.4662620 551.22 8.212435 10.044170 b
## Northeast 9.448888 0.4707121 559.55 8.524310 10.373467 b
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05posthoc<-emmeans(mod, ~Season)
CLD(posthoc, Letters=letters)## Season emmean SE df lower.CL upper.CL .group
## wet season 7.771551 0.3348801 206.76 7.111333 8.431768 a
## dry season 12.181502 0.3071824 206.48 11.575885 12.787118 b
##
## Results are averaged over the levels of: Location
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05
Figure DLI model output
PAR and DLI at 3 depths
With Depth as a factor (<1m, 2m, 8m), we can produce a continuous graph of light (PAR or DLI) at 2m and use model intercepts to predict the light at the shallow and deep zones (ca. <1m and 8m, respectfully). Here, the model is \(lm(log.PAR\) ~ \(Depth.factor)\).
# Plot it!!
# Reef 44
par(mfrow=c(2,2), mar=c(4,5,2,2))
plot(y=All.PAR.df$Rf44.shall, x=All.PAR.df$timestamp, type="l", col="palegreen3",
cex.main=0.7, ylab=(expression(paste("PAR"~(mu*mol~photons~m^-2~s^-1), sep="")))
, xlab="", main="Northwest", cex.main=1, ylim=c(0, 2400))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf44.mid, x=All.PAR.df$timestamp, type="l", col="gray"))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf44.deep, x=All.PAR.df$timestamp, type="l", col="palegreen4"))
legend("topright", lty=1, col=c("palegreen3", "gray", "palegreen4"), legend=c("<1m", "2m", "8m"),
lwd=2, bty="n", inset=c(0, -0.01), seg.len=1, cex=0.9, x.intersp=0.5, y.intersp=0.8)
# Plot it!!
# Reef 42
plot(y=All.PAR.df$Rf42.shall, x=All.PAR.df$timestamp, type="l", col="mediumturquoise",
cex.main=0.7, ylab=(expression(paste("PAR"~(mu*mol~photons~m^-2~s^-1), sep="")))
, xlab="", main="Northeast", cex.main=1, ylim=c(0, 2400))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf42.mid, x=All.PAR.df$timestamp, type="l", col="gray"))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf42.deep, x=All.PAR.df$timestamp, type="l", col="dodgerblue2"))
legend("topright", lty=1, col=c("mediumturquoise", "gray", "dodgerblue2"), legend=c("<1m", "2m", "8m"),
lwd=2, bty="n", inset=c(0, -0.01), seg.len=1, cex=0.9, x.intersp=0.5, y.intersp=0.8)
# Plot it!!
# Reef 10
plot(y=All.PAR.df$Rf10.shall, x=All.PAR.df$timestamp, type="l", col="orangered",
cex.main=0.7, ylab=(expression(paste("PAR"~(mu*mol~photons~m^-2~s^-1), sep="")))
, xlab="Dates", main="Southwest", cex.main=1, ylim=c(0, 2400))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf10.mid, x=All.PAR.df$timestamp, type="l", col="gray"))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$Rf10.deep, x=All.PAR.df$timestamp, type="l", col="lightsalmon3"))
legend("topright", lty=1, col=c("orangered", "gray", "lightsalmon3"), legend=c("<1m", "2m", "8m"),
lwd=2, bty="n", inset=c(0, -0.01), seg.len=1, cex=0.9, x.intersp=0.5, y.intersp=0.8)
# Plot it!
# HIMB
plot(y=All.PAR.df$HIMB.shall, x=All.PAR.df$timestamp, type="l", col="mediumorchid2",
cex.main=0.7, ylab=(expression(paste("PAR"~(mu*mol~photons~m^-2~s^-1), sep="")))
, xlab="Dates", main="Southeast", cex.main=1, ylim=c(0, 2400))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$HIMB.mid, x=All.PAR.df$timestamp, type="l", col="gray"))
par(new=T)
with(All.PAR.df, lines(y=All.PAR.df$HIMB.deep, x=All.PAR.df$timestamp, type="l", col="plum4"))
legend("topright", lty=1, col=c("mediumorchid2", "gray", "plum4"), legend=c("<1m", "2m", "8m"),
lwd=2, bty="n", inset=c(0, -0.01), seg.len=1, cex=0.9, x.intersp=0.5, y.intersp=0.8)
Figure PAR at the four reef areas from summer to winter 2016 at three depths (<1m, 2m, 8m). Data at 2m was recorded from loggers; data at <1m and 8m were modeled using light attenuation calculations
dev.copy(pdf, "figures/environmental/PAR.all depths.pdf", height=7, width=8)
dev.off()#######
####### calculate DLI and compare
#######
df<-All.PAR.df
df$timestamp<-strptime(df$timestamp, format="%Y-%m-%d %H:%M:%S")
df$timestamp<-as.Date(df$timestamp, format="%Y-%m-%d") # change to DATE ONLY format to calulate range, means
df[is.na(df)] <- 0
df.split <- split(df, f=df$timestamp
< as.Date("2016-06-10", format="%Y-%m-%d")) # split df by date
df.dli<-aggregate(data.frame(Rf42.shall.DLI=df.split[[1]]$Rf42.shall*0.0864,
Rf42.mid.DLI=df.split[[1]]$Rf42.mid*0.0864,
Rf42.deep.DLI=df.split[[1]]$Rf42.deep*0.0864,
Rf44.shall.DLI=df.split[[1]]$Rf44.shall*0.0864,
Rf44.mid.DLI=df.split[[1]]$Rf44.mid*0.0864,
Rf44.deep.DLI=df.split[[1]]$Rf44.deep*0.0864,
HIMB.shall.DLI=df.split[[1]]$HIMB.shall*0.0864,
HIMB.mid.DLI=df.split[[1]]$HIMB.mid*0.0864,
HIMB.deep.DLI=df.split[[1]]$HIMB.deep*0.0864,
Rf10.shall.DLI=df.split[[1]]$Rf10.shall*0.0864,
Rf10.mid.DLI=df.split[[1]]$Rf10.mid*0.0864,
Rf10.deep.DLI=df.split[[1]]$Rf10.deep*0.0864),
by=list(Date=df.split[[1]]$timestamp), FUN=mean)
# make a date sequence for entire study
all.date.time<-as.data.frame(seq(
from=as.POSIXct("2016-06-10", tz="HST"),
to=as.POSIXct("2017-01-12", tz="HST"),
by="1 d"))
colnames(all.date.time)[1]<-"Date"
all.date.time$Date<-strptime(all.date.time$Date, format="%Y-%m-%d")
all.date.time$Date<-as.Date(all.date.time$Date, format="%Y-%m-%d")
df.dli<-merge(all.date.time, df.dli, by="Date", all.x=T)
df.dli[df.dli<=0] <- NA
#####################
# Plot it!!
# Reef 44
par(mfrow=c(2,2), mar=c(4,5,2,2))
plot(y=df.dli$Rf44.shall, x=df.dli$Date, type="h", col="palegreen3",
cex.main=0.7, ylab=(expression(paste("DLI"~(mol~photons~m^-2~d^-1), sep="")))
, xlab="", main="Northwest", cex.main=1, ylim=c(0, 50))
par(new=T)
with(df.dli, lines(y=df.dli$Rf44.mid, x=df.dli$Date, type="h", col="gray"))
par(new=T)
with(df.dli, lines(y=df.dli$Rf44.deep, x=df.dli$Date, type="h", col="palegreen4"))
legend("topright", lty=1, col=c("palegreen3", "gray", "palegreen4"), legend=c("<1m", "2m", "8m"),
lwd=2, bty="n", inset=c(0.01, 0), seg.len=1, cex=0.9, x.intersp=0.2, y.intersp=1)
# Plot it!!
# Reef 42
plot(y=df.dli$Rf42.shall, x=df.dli$Date, type="h", col="mediumturquoise",
cex.main=0.7, ylab=(expression(paste("DLI"~(mol~photons~m^-2~d^-1), sep="")))
, xlab="", main="Northeast", cex.main=1, ylim=c(0, 50))
par(new=T)
with(df.dli, lines(y=df.dli$Rf42.mid, x=df.dli$Date, type="h", col="gray"))
par(new=T)
with(df.dli, lines(y=df.dli$Rf42.deep, x=df.dli$Date, type="h", col="dodgerblue2"))
legend("topright", lty=1, col=c("mediumturquoise", "gray", "dodgerblue2"), legend=c("<1m", "2m", "8m"),
lwd=2, bty="n", inset=c(0.01, 0), seg.len=1, cex=0.9, x.intersp=0.2, y.intersp=1)
# Plot it!!
# Reef 10
plot(y=df.dli$Rf10.shall, x=df.dli$Date, type="h", col="orangered",
cex.main=0.7, ylab=(expression(paste("DLI"~(mol~photons~m^-2~d^-1), sep="")))
, xlab="Dates", main="Southwest", cex.main=1, ylim=c(0, 50))
par(new=T)
with(df.dli, lines(y=df.dli$Rf10.mid, x=df.dli$Date, type="h", col="gray"))
par(new=T)
with(df.dli, lines(y=df.dli$Rf10.deep, x=df.dli$Date, type="h", col="lightsalmon3"))
legend("topright", lty=1, col=c("orangered", "gray", "lightsalmon3"), legend=c("<1m", "2m", "8m"),
lwd=2, bty="n", inset=c(0.01, 0), seg.len=1, cex=0.9, x.intersp=0.2, y.intersp=1)
# Plot it!
# HIMB
plot(y=df.dli$HIMB.shall, x=df.dli$Date, type="h", col="mediumorchid2",
cex.main=0.7, ylab=(expression(paste("DLI"~(mol~photons~m^-2~d^-1), sep="")))
, xlab="Dates", main="Southeast", cex.main=1, ylim=c(0, 50))
par(new=T)
with(df.dli, lines(y=df.dli$HIMB.mid, x=df.dli$Date, type="h", col="gray"))
par(new=T)
with(df.dli, lines(y=df.dli$HIMB.deep, x=df.dli$Date, type="h", col="plum4"))
legend("topright", lty=1, col=c("mediumorchid2", "gray", "plum4"), legend=c("<1m", "2m", "8m"),
lwd=2, bty="n", inset=c(0.01, 0), seg.len=1, cex=0.9, x.intersp=0.2, y.intersp=1)
Figure DLI at the four reef areas from summer to winter 2016 at three depths (<1m, 2m, 8m). Data at 2m was recorded from loggers; data at <1m and 8m were modeled using light attenuation calculations
dev.copy(pdf, "figures/environmental/DLIcalc.all depths.pdf", height=7, width=8)
dev.off()DLI.summary<-read.csv("data/environmental/temp and light/Jun_DecPAR/DLI.3depth.summary.csv")
DLI.summary$Location<-factor(DLI.summary$Location, levels=c("NW", "NE", "SW", "SE"))
DLI.summary$Depth.zone<-factor(DLI.summary$Depth.zone, levels=c("shallow", "mid", "deep"))
#####
# figure formatting script
Fig.formatting<-(theme_classic()) +
theme(text=element_text(size=10),
axis.line=element_blank(),
legend.justification=c(1,1), legend.position=c(1,1),
legend.background = element_rect("transparent", colour=NA),
legend.text=element_text(size=10),
legend.title = element_blank(),
panel.border = element_rect(fill=NA, colour = "black", size=0.8),
aspect.ratio=1,
axis.ticks = element_line(colour = 'black', size = 0.4),
axis.ticks.length=unit(0.25, "cm"),
axis.text.y=element_text(
margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=10),
axis.text.x=element_text(
margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=10)) +
theme(legend.spacing.x = unit(0.2, 'cm'),
legend.key.size = unit(0.4, "cm"))
pd <- position_dodge(0.7) #offset for error bars and columns
#####
Location.label<-c("NW", "NE", "SW", "SE")
Depth.label<-c("<1m", "2m", "8m")
plot_colors2<-c("lightskyblue", "dodgerblue", "deepskyblue4")
plot_stat<-c("slategray", "dodgerblue2", "steelblue1", "slategray", "dodgerblue2", "steelblue1",
"slategray", "dodgerblue2", "steelblue1","slategray", "dodgerblue2", "steelblue1")
# DLI by site and depth
Fig.DLI<-ggplot(DLI.summary, aes(x=Location, y=DLI.mean, fill=Depth.zone)) +
geom_bar(stat="identity", col=plot_stat, position = pd, width=0.7, size=0.3) +
geom_errorbar(aes(ymin=DLI.mean-DLI.se, ymax=DLI.mean+DLI.se),
size=0.4, width=0, position=pd, col=plot_stat) +
xlab("Locations") +
ylab(expression(paste("DLI", ~(mol~photons~m^-2~d^-1), sep=""))) +
scale_x_discrete(labels=Location.label) +
scale_y_continuous(expand=c(0,0), limits=c(0, 35)) +
scale_fill_manual(values=alpha(c(plot_colors2), 0.6),
labels=Depth.label)+ Fig.formatting; Fig.DLI
Figure Daily Light Integral (DLI) for four reef locations at three depths. DLI estimated from light attenuation coefficents using depth as a fixed effect against baseline loggers at 2 m depth
ggsave(file="figures/environmental/DLI.bar.pdf", height=4, width=4)Attenuation and Light-at-Depth
A periodic deployment of loggers (October and November) at ca. <1m and ca. 7m were used to model the relationship between 2m light and light at other depths.
Using the logger data, the difference in light at 1 and 7m was calculated relative to 2m (i.e, PAR baseline (2m) - PAR other), then the difference in depth at 1 and 7m was calculated realtive to logger depth i.e, depth logger other - depth logger (2m baseline). A no-intercept model was used to fit difference in log(DLI) and depth (i.e, log(delta-DLI) and delta-depth) with depth as a continous/numeric variable to determine kd: \(lm(delta.log.DLI\)~\(delta.Depth + 0)\).
Site-specific kd values were used to generate a continuous variable of “light at depth”. At each time period (summer or winter) the daily light integral was averaged over 3 months: summer 2016 (June, July, August) and winter 2016 (November, December, January). Discrete monthly seasonal-mean DLI means was used to generate an estimate of light at depth for corals at each site, as this integrates the site and seasoanl variance.
Light at Depth was calculated following Beer-Lambert: \(Ez{_d} = Ez_{_0}e^{(-k_d *d)}\), where \(d\) is depth, \(k_d\) is the attenuation coefficient for each site, \(Ez{_d}\) is light at depth \(d\) and \(Ez{_0}\) is measured light. Calculations were relative to light at 2m depth (i.e., \(Ez_{_0}\) was the light at 2m depth) and \(d\) was the difference between coral depth relative to depth of the logger (i.e., 2m). light reaching depth d and Ez(0) is the incident light at the surface.
########################## # LIGHT AT DEPTH
#############
######
# using new depth and the site-specific kd (light attenuation) determine approximate "light at depth" for each colony sample.
######
### attach necessary files
logger.depths<-read.csv("data/environmental/temp and light/light.logger.depths.csv") # depths for loggers
kd.all<-read.csv("data/environmental/kd.all.csv") # kds for each site
# Seasonal DLI used for "period of collection" light levels
month.DLI<-read.csv("data/environmental/temp and light/Jun_DecPAR/All.DLI_long.csv")
# corals collected in June-July-August use summer time DLI for these months as indicator of average DLI
summer.DLI<-month.DLI[(month.DLI$Month=="June" | month.DLI$Month=="July" | month.DLI$Month=="August"),]
winter.DLI<-month.DLI[(month.DLI$Month=="November" | month.DLI$Month=="December" |month.DLI$Month=="January"),]
# summer mean and SE dataframe
sum.mean<-aggregate(DLI~Site, summer.DLI, mean)
sum.SE<-aggregate(DLI~Site, summer.DLI, std.error)
sum.light.df<-cbind(sum.mean, sum.SE[2])
colnames(sum.light.df)=c("Site", "mean.DLI", "se.DLI")
sum.light.df$Season<-as.factor("summer")
# winter mean and SE dataframe
wint.mean<-aggregate(DLI~Site, winter.DLI, mean)
wint.SE<-aggregate(DLI~Site, winter.DLI, std.error)
wint.light.df<-cbind(wint.mean, wint.SE[2])
colnames(wint.light.df)=c("Site", "mean.DLI", "se.DLI")
wint.light.df$Season<-as.factor("winter")
season.DLI<-rbind(sum.light.df[,c(4,1,2:3)], wint.light.df[,c(4,1,2:3)]) # compiled means for DLI at ~2m
season.DLI$Location<- as.factor(c("SE", "SW", "NE", "NW", "SE", "SW", "NE", "NW")); season.DLI<-season.DLI[,c(1,5,2:4)]
write.csv(season.DLI, "data/environmental/season.DLI.csv")
#######################################
### make new dataframe for calculations
df.light<- data[, c("Season", "Reef", "Location", "newDepth")]
# make a column of "depth differences" relative to where ~2m logger was deployed
df.light$depth.diff<-ifelse(df.light$Reef=="F1-R46", df.light$newDepth-1.8288,
ifelse(df.light$Reef=="F8-R10", df.light$newDepth-1.8288,
ifelse(df.light$Reef=="HIMB", df.light$newDepth-1.8288,
df.light$newDepth-2.4384))) # last statement for remaining site (Rf42)
# make a column for sample-specific light at depth (estimate) based on kd
# follow: #with 2m as baseline# E(depth) = E(2m)*exp(-k_d*(depth - 2m))
# so that: light at depth = DLI at mid.depth * exp (-kd *(delta shallow-deep in m))
df.light$Light<-ifelse(df.light$Reef=="HIMB" & df.light$Season=="summer",
season.DLI$mean.DLI[1]*exp(-kd.all$kd[1]*df.light$depth.diff), # summer HIMB
ifelse(df.light$Reef=="F8-R10" & df.light$Season=="summer",
season.DLI$mean.DLI[2]*exp(-kd.all$kd[2]*df.light$depth.diff), # summer R10
ifelse(df.light$Reef=="R42" & df.light$Season=="summer",
season.DLI$mean.DLI[3]*exp(-kd.all$kd[3]*df.light$depth.diff), # summer R42
ifelse(df.light$Reef=="F1-R46" & df.light$Season=="summer",
season.DLI$mean.DLI[4]*exp(-kd.all$kd[4]*df.light$depth.diff), # summer R46
ifelse(df.light$Reef=="HIMB" & df.light$Season=="winter",
season.DLI$mean.DLI[5]*exp(-kd.all$kd[1]*df.light$depth.diff), # winter HIMB
ifelse(df.light$Reef=="F8-R10" & df.light$Season=="winter",
season.DLI$mean.DLI[6]*exp(-kd.all$kd[2]*df.light$depth.diff), # winter R10
ifelse(df.light$Reef=="R42" & df.light$Season=="winter",
season.DLI$mean.DLI[7]*exp(-kd.all$kd[3]*df.light$depth.diff), # winter R42
season.DLI$mean.DLI[8]*exp(-kd.all$kd[4]*df.light$depth.diff)) # winter R46
))))))
###### plot of light x depth by season
df.light$Reef <- factor(df.light$Reef, levels = c("F1-R46", "R42", "F8-R10", "HIMB"))
df.light$Location <- factor(df.light$Location, levels=c("NW", "NE", "SW", "SE"))
plot.by.sites=ggplot(df.light, aes(Light)) + geom_density(aes(fill=Location), alpha=0.3, position = 'stack') + scale_x_continuous(limits=c(0, 40)) + ggtitle("Light at Depth by Seasons")
plot.by.site=ggplot(df.light, aes(Light)) + geom_density(aes(fill=Season), alpha=0.3, position = 'stack') +
scale_x_continuous(limits=c(0, 40)) + ggtitle("Light at Depth by Seasons") + facet_wrap(~Location, scales="free") + scale_fill_manual(values=c("darkorange1", "dodgerblue1"))
######
# can loop as a list
p=ggplot(df.light, aes(Light, fill=Season)) + geom_density(alpha=0.3, position = 'stack') + scale_x_continuous(limits=c(0, 40)) + ggtitle("Light at Depth by Seasons")
plots=dlply(df.light, .(Location), function(x) p %+% x + facet_wrap(~Location))
###### plot of light x depth by season with exponential curve fitting
Sum<-df.light[(df.light$Season=="summer"),]
Win<-df.light[(df.light$Season=="winter"),]
plot(Light~newDepth, Sum, col="coral", pch=16, xlab="Depth (m)", ylab="DLI at coral",
ylim=c(0, 30), xlim=c(0, 10))
summod<-lm(log(Light)~newDepth, Sum)
curve(exp(coef(summod)["(Intercept)"]+coef(summod)["newDepth"]*x), add=TRUE, col="coral", lwd=2)
par(new=T)
plot(Light~newDepth, Win, col="lightblue2", pch=16, xaxt="n", yaxt="n",
xlab="", ylab="", ylim=c(0, 30), xlim=c(0, 10))
wintmod<-lm(log(Light)~newDepth, Win)
curve(exp(coef(wintmod)["(Intercept)"]+coef(wintmod)["newDepth"]*x), add=TRUE, col="lightblue2", lwd=2)
legend("topright", legend=c("summer", "winter"), col=c("coral", "lightblue2"), lty=1, lwd=1, pch=16, bty="n")
Figure Relationship between mean seasonal light availability (mean daily light integral) at depth of coral fragment collection
Dissolved nutrients
Dissolved inorganic nutrients in seawater were quantified on water samples (ca. 100 ml) taken from each site and filtered through a GF/F filter and stored in acid-washed bottles and frozen at -20 °C until analyzes. Dissolved inorganic nutrients—ammonium (NH4+), nitrate + nitrite (NO3- + NO2-), phosphate (PO43-), and silicate (Si(OH)4)—were measured at University of Hawai‘i at Mānoa SOEST Laboratory for Analytical Biogeochemistry using a Seal Analytical AA3 HR nutrient autoanalyzer and expressed as μmol l-1.
Silicate showed no effects (p > 0.323). Phosphate increased in the winter (p=0.046) most notably at northern sites. N+N was higher in northern sites relative to southern sites and increased during the winter at all sites compared to the summer. Ammonium increased in the winter as well (p <0.001).
#########################
#########################
# nutrients
nutr<-read.csv("data/environmental/PanKBay_nutrients.csv")
#models
mod<-lm(silicate~Location+Date, data=nutr); Anova(mod, type=2)
mod<-lm(phosphate~Location+Date, data=nutr); Anova(mod, type=2)
mod<-lm(N.N~Location+Date, data=nutr); Anova(mod, type=2)
mod<-lm(ammonium~Location+Date, data=nutr); Anova(mod, type=2)
nutr$Date<-mdy(nutr$Date) # corrects date format
dates<-cbind("10-Aug '16", "20-Dec '16")
RfHIMB<-nutr[(nutr$Reef=="HIMB"), ]
Rf42<-nutr[(nutr$Reef=="R42"), ]
Rf46<-nutr[(nutr$Reef=="F1-46"), ]
RF10<-nutr[(nutr$Reef=="F8-R10"), ]
Reefs<-c("Reef 46", "Reef 42", "Reef 10", "HIMB")
Locations<-c("NW", "NE", "SW", "SE")
plot_colors<-c("mediumseagreen", "dodgerblue", "salmon", "orchid")
###############
# Phosphate
# par(mfrow=c(1,1), mar=c(2,4,1,1), mgp=c(2,0.5,0))
# figure layout
layout(matrix(c(1,2,3,4), nrow=1, byrow=TRUE))
par(mar=c(5,4.5,7,1.5))
###############
# Silicate
plot(y=Rf46$silicate, x=Rf46$Date, xaxt="n", type="o", xlab=NA, ylab=expression(paste("silicate"~(mu*mol~L^-1), sep="")), ylim=c(0, 15), pch=19, lty=2, cex=1, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis(side=1, at=Rf46$Date, labels=dates, cex.axis=0.8) # plots Reef 46 / NW
#plots Reef 42/ NE
with(Rf46, lines(Rf42$silicate, x=Rf42$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[2]))
#plots Reef 10/ SW
with(Rf46, lines(RF10$silicate, x=RF10$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[3]))
#plots HIMB/ SE
with(RfHIMB, lines(RfHIMB$silicate, x=RfHIMB$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[4]))
###############
# Phosphate
plot(y=Rf46$phosphate, x=Rf46$Date, xaxt="n", type="o", xlab=NA, ylab=expression(paste("phosphate"~(mu*mol~L^-1), sep="")), ylim=c(0, 0.25), pch=19, lty=2, cex=1, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis(side=1, at=Rf46$Date, labels=dates, cex.axis=0.8) # plots Rf46 /NW
#plots Reef 42 / NE
with(Rf46, lines(Rf42$phosphate, x=Rf42$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[2]))
#plots Reef 10 /SW
with(Rf46, lines(RF10$phosphate, x=RF10$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[3]))
#plots HIMB / SE
with(RfHIMB, lines(RfHIMB$phosphate, x=RfHIMB$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[4]))
###############
# N+N
plot(y=Rf46$N.N, x=Rf46$Date, xaxt="n", type="o", xlab=NA, ylab=expression(paste("nitrate+nitrite"~(mu*mol~L^-1), sep="")), ylim=c(0, 1.5), pch=19, lty=2, cex=1, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis(side=1, at=Rf46$Date, labels=dates, cex.axis=0.8) # plots Reef 46/ NW
#plots Reef 42/ NE
with(Rf46, lines(Rf42$N.N, x=Rf42$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[2]))
#plots Reef 10 / SE
with(Rf46, lines(RF10$N.N, x=RF10$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[3]))
#plots HIMB / SW
with(RfHIMB, lines(RfHIMB$N.N, x=RfHIMB$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[4]))
###############
# Ammonium
plot(y=Rf46$ammonium, x=Rf46$Date, xaxt="n", type="o", xlab=NA, ylab=expression(paste("ammonium"~(mu*mol~L^-1), sep="")), ylim=c(0, 1.5), pch=19, lty=2, cex=1, lwd=1, col=plot_colors[1], cex.axis=0.8)
axis(side=1, at=Rf46$Date, labels=dates, cex.axis=0.8) # plots Reef 46/ NW
#plots Reef 42 /NE
with(Rf46, lines(Rf42$ammonium, x=Rf42$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[2]))
#plots Reef 10/ SW
with(Rf46, lines(RF10$ammonium, x=RF10$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[3]))
#plots HIMB /SE
with(RfHIMB, lines(RfHIMB$ammonium, x=RfHIMB$Date, type="o", pch=19, lty=2, cex=1, lwd=1, col=plot_colors[4]))
legend("topleft", legend=Locations, col=plot_colors[1:4], lwd=1, pch=19, lty=2, cex=0.8, bty="n", x.intersp=0.2, y.intersp=1.2, inset=c(0.01, 0), seg.len=1.5)
dev.copy(pdf, "figures/environmental/all.nutrients.pdf", encod="MacRoman", height=3, width=7)
dev.off()Suspended particulates
Gravimetric SPM
Total (organic + inorganic) suspended particulate matter (SPM) was quantified by collecting and filtering 0.5 L of seawater from each location during each season. Seawater was sieved at 243 μm to remove large debris and filtered onto a pre-combusted (4h, 450°C) GF/F filter (0.7 μm). Salts were removed by light rinsing with ddH2O. Filters were dried at 60°C for 24h and re-weighed to closest 0.001 g to obtain total SPM. Masses were normalized to a liter of water.
###### SPM
spm<-read.csv("data/environmental/PanKBay_SPM.csv")
spm$SPM..g.l<-spm$SPM..g.l*1000 # change to mg
spm$Date<-factor(spm$Date, levels=c("8/10/16", "12/20/16"))
colnames(spm)[4]="SPM.mg"
spm$Location<-ordered(spm$Location, c("NW", "NE", "SW", "SE"))
mod<-lm(SPM.mg~Location+Date, data=spm); Anova(mod, type=2)## Anova Table (Type II tests)
##
## Response: SPM.mg
## Sum Sq Df F value Pr(>F)
## Location 8.8382 3 10.3062 0.0003025 ***
## Date 1.4114 1 4.9373 0.0386252 *
## Residuals 5.4312 19
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1posthoc<-emmeans(mod, ~Location)
CLD(posthoc, Letters=letters)## Location emmean SE df lower.CL upper.CL .group
## NE 0.7966668 0.2182711 19 0.3398203 1.253513 a
## NW 1.6450002 0.2182711 19 1.1881536 2.101847 ab
## SE 2.2666668 0.2182711 19 1.8098203 2.723513 b
## SW 2.2850002 0.2182711 19 1.8281536 2.741847 b
##
## Results are averaged over the levels of: Date
## Confidence level used: 0.95
## P value adjustment: tukey method for comparing a family of 4 estimates
## significance level used: alpha = 0.05plot(allEffects(mod), par.strip.text=list(cex=0.7))
rm(mean) # this will remove 'mean' from library and let the function below proceed.
mean<-aggregate(SPM.mg~Date + Location+Reef, spm, mean)
se<-aggregate(SPM.mg~Date+Location+Reef, spm, std.error)
spm.df<-cbind(mean, se[c(0,4)])
colnames(spm.df)=c("Date", "Location", "Reef", "SPM", "se")
spm.df$Date<-ordered(spm.df$Date, c("8/10/16", "12/20/16"))
spm.df$Reef<-ordered(spm.df$Reef, c("F1-46", "R42", "F8-10", "HIMB"))
spm.df$Location<-ordered(spm.df$Location, c("NW", "NE", "SW", "SE"))
Reefs<-c("Reef 46", "Reef 42", "Reef 10", "HIMB")
Locations<-c("NW", "NE", "SW", "SE")
Date<-c("10-Aug '16", "20-Dec '16")
plot_colors2<-c("mediumturquoise", "skyblue3", "lightcoral", "plum")
#####
# figure formatting script
Fig.formatting<-(theme_classic()) +
theme(text=element_text(size=10),
axis.line=element_blank(),
legend.justification=c(1,1), legend.position=c(1,1),
legend.background = element_rect("transparent", colour=NA),
legend.text=element_text(size=10),
legend.title = element_blank(),
panel.border = element_rect(fill=NA, colour = "black", size=0.8),
aspect.ratio=1,
axis.ticks = element_line(colour = 'black', size = 0.4),
axis.ticks.length=unit(0.25, "cm"),
axis.text.y=element_text(
margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=10),
axis.text.x=element_text(
margin=unit(c(0.5, 0.5, 0.5, 0.5), "cm"), colour="black", size=10)) +
theme(legend.spacing.x = unit(0.2, 'cm'),
legend.key.size = unit(0.4, "cm"))
pd <- position_dodge(0.7) #offset for error bars and columns
#####
# SPM (mg/l)
Fig.SPM<-ggplot(spm.df, aes(x=Date, y=SPM, fill=Location)) +
geom_bar(colour=c("plum","lightcoral", "skyblue3","mediumturquoise","plum","lightcoral", "skyblue3", "mediumturquoise"),
stat="identity", position = pd, width=0.7, size=0.3) +
geom_errorbar(aes(ymin=SPM-se, ymax=SPM+se),
colour=c("plum","lightcoral", "skyblue3","mediumturquoise","plum","lightcoral", "skyblue3","mediumturquoise"),
size=0.4, width=0, position=pd) +
xlab("Sampling points") +
ylab(expression(paste("SPM", ~(mg~L^-1), sep=""))) +
scale_x_discrete(labels=Date) +
scale_y_continuous(expand=c(0,0), limits=c(0, 5)) +
scale_fill_manual(values=alpha(c(plot_colors2), 0.5),
labels=Locations)+ Fig.formatting; Fig.SPM
ggsave(file="figures/environmental/SPM.pdf", height=4, width=4)Isotopes of SPM/plankton
Plankton tows and seawater collections to parameterized particulates and plankters as potential heterotrophic end members available to reef corals. Sampling was done at 4 sites where corals were collected [North: (Reef 42, fringe-Reef 46), and South: (HIMB, fringe-Reef 10)], as well as two reefs in central region of the bay where corals were not collected (Central: Reef 25, fringe-Reef 25) to increase spatial resolution of suspended particulate isotope sample sizes. Using size-fractioned materials collected in seawater and plankton tows, we examined the spatial (among reef sites) and temporal patterns (among seasons) in stable isotope values of size-fractioned end members. Carbon and nitrogen isotope values among the 6 reef sites were not significantly different (p > 0.140), season had marginal effects (p > 0.049), whereas fraction influenced isotope values significantly (p< 0.001). Therefore, isotope values were pooled among reefs and seasons to best interpret size-fraction isotope values. Samples were not acidified prior to analysis as to not affect nitrogen isotope values.
SWiso<-read.csv("data/isotopes_SW_all times.csv")
mod<-(lm(d13C~Location+Season+SW.fraction..um, data=SWiso)); Anova(mod, type=2) # no location effect## Anova Table (Type II tests)
##
## Response: d13C
## Sum Sq Df F value Pr(>F)
## Location 17.926 5 1.3418 0.2626013
## Season 7.921 1 2.9645 0.0914188 .
## SW.fraction..um 76.419 4 7.1504 0.0001286 ***
## Residuals 130.920 49
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1mod<-(lm(d15N~Location+Season+SW.fraction..um, data=SWiso)); Anova(mod, type=2) # no location effect## Anova Table (Type II tests)
##
## Response: d15N
## Sum Sq Df F value Pr(>F)
## Location 2.3255 5 1.7291 0.14557
## Season 1.0935 1 4.0653 0.04927 *
## SW.fraction..um 30.3773 4 28.2335 3.454e-12 ***
## Residuals 13.1802 49
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1posthoc<-emmeans(mod, ~Season)
CLD(posthoc, Letters=letters)## Season emmean SE df lower.CL upper.CL .group
## summer 6.25 0.09468949 49 6.059714 6.440286 a
## winter 6.52 0.09468949 49 6.329714 6.710286 b
##
## Results are averaged over the levels of: Location, SW.fraction..um
## Confidence level used: 0.95
## significance level used: alpha = 0.05SWiso$Reef<-factor(SWiso$Reef, levels=c("F1-46", "R42", "F2-R25", "R25", "F8-R10", "HIMB"))
SWiso$Location<-factor(SWiso$Location, levels=c("NW", "NE", "CW", "CE", "SW", "SE"))
SWiso$SW.fraction..um<-factor(SWiso$SW.fraction..um, levels=c(">243", "<243", "100-243", "10-100", "0-10"))
winter.data<-SWiso[(SWiso$Season=="winter"),]
summer.data<-SWiso[(SWiso$Season=="summer"),]These box plots show the seasonal isotope values for carbon and nitrogen according to size fractions.
# plots of SW isotopes by Season-- first showing by size, then by site
op<-par(mfrow = c(2,2), mar=c(5,5,2,1))
######### Sizes
# Summer size d13C
plot(d13C~SW.fraction..um, data=summer.data, ylim=c(-25,-10),
main=expression(paste("summer"~ delta^{13}, "C")), ylab=expression(paste(delta^{13}, "C (‰, V-PDB)")),
col="paleturquoise3", cex.axis=0.7, cex.main=1, cex.lab= 0.8, xlab=expression(paste("Size Fraction (", mu,m,")")))
# Winter size d13C
plot(d13C~SW.fraction..um, data=winter.data, ylim=c(-25,-10),
main=expression(paste("winter"~ delta^{13}, "C")), ylab=expression(paste(delta^{13}, "C (‰, V-PDB)")),
col="paleturquoise3", cex.axis=0.7, cex.main=1, cex.lab= 0.8, xlab=expression(paste("Size Fraction (", mu,m,")")))
# Summer size d15N
plot(d15N~SW.fraction..um, data=summer.data, ylim=c(3,10),
main=expression(paste("summer"~ delta^{15}, "N")), col="plum", cex.axis=0.7, cex.main=1, cex.lab= 0.8,
xlab=expression(paste("Size Fraction (", mu,m,")")),
ylab=expression(paste(delta^{15}, "N (‰, air)")))
# Winter size d15N
plot(d15N~SW.fraction..um, data=winter.data, ylim=c(3,10),
main=expression(paste("winter"~ delta^{15}, "N")), ylab=expression(paste(delta^{15}, "N (‰, air)")),
col="plum", cex.axis=0.7, cex.main=1, cex.lab= 0.8, xlab=expression(paste("Size Fraction (", mu,m,")")))
Figure Size-fractioned sample isotope values stratified by seasons (summer and winter)
These figures show the same size-fractioned isotope data as above, but stratified by location and pooled across seasons.
######## Sites
op<-par(mfrow = c(2,2), mar=c(5,5,2,1))
# Summer site d13C
plot(d13C~Location, data=summer.data, ylim=c(-25,-10),
main=expression(paste("summer"~ delta^{13}, "C")), ylab=expression(paste(delta^{13}, "C (‰, V-PDB)")),
col="lightskyblue", cex.axis=0.7, cex.main=1, cex.lab= 0.8, xlab="Reef Sites")
# Winter site d13C
plot(d13C~Location, data=winter.data, ylim=c(-25,-10),
main=expression(paste("winter"~ delta^{13}, "C")), ylab=expression(paste(delta^{13}, "C (‰, V-PDB)")),
col="lightskyblue", cex.axis=0.7, cex.main=1, cex.lab= 0.8, xlab="Reef Sites")
# Summer site d15N
plot(d15N~Location, data=summer.data, ylim=c(3,10),
main=expression(paste("summer"~ delta^{15}, "N")), col="salmon", cex.axis=0.7, cex.main=1, cex.lab= 0.8,
ylab=expression(paste(delta^{15}, "N (‰, air)")), xlab="Reef Sites")
# Winter site d15N
plot(d15N~Location, data=winter.data, ylim=c(3,10),
main=expression(paste("winter"~ delta^{15}, "N")), ylab=expression(paste(delta^{15}, "N (‰, air)")),
col="salmon", cex.axis=0.7, cex.main=1, xlab="Reef Sites")
Figure Isotope values in seawater particles stratified by sites (northwest to southeast)
Pooling the isotope data across sites and seasons (where effects were absent or negligible), this planktonic foodweb for carbon and nitrogen illustrates differences in the food sources and their isotopic composition.
mod<-lm(d13C~SW.fraction..um, data=SWiso); Anova(mod, type=2)
posthoc<-emmeans(mod, ~SW.fraction..um)
CLD(posthoc, Letters=letters)
mod<-lm(d15N~SW.fraction..um, data=SWiso); Anova(mod, type=2)
posthoc<-emmeans(mod, ~SW.fraction..um)
CLD(posthoc, Letters=letters)
rm(mean)
####
#### making scatter for d15N and d13C, pooled across seasons and sites
mix.N.mean<-aggregate(d15N~SW.fraction..um, data=SWiso, mean)
mix.N.se<-aggregate(d15N~SW.fraction..um, data=SWiso, std.error)
mix.C.mean<-aggregate(d13C~SW.fraction..um, data=SWiso, mean)
mix.C.se<-aggregate(d13C~SW.fraction..um, data=SWiso, std.error)
mix.data<-cbind(mix.N.mean, mix.C.mean[c(2,0)], mix.N.se[c(2,0)], mix.C.se[c(2,0)]); colnames(mix.data)=c("fraction", "d15N", "d13C", "d15N.se", "d13C.se")
colors=c("#FF6A6A", "#00B2EE", "#FFB90F", "#3CB371", "#8B7500")
op<-par(mfrow = c(1,1), mar=c(5,4,1,5),xpd=TRUE, pty="sq")
size.labels=c(expression(paste("> 243"~mu,m)), expression(paste("< 243"~mu,m)), expression(paste("100-243"~mu,m)), expression(paste("10-100"~mu,m)), expression(paste("< 10"~mu,m)))
#### make the plot
plot(d15N~d13C, data=mix.data, type="n", ylim=c(4.5,8), xlim=c(-23,-17), tck=-0.03, cex.axis=0.7, cex.lab=0.7,
ylab=expression(paste(delta^{15}, "N (‰, air)")), xlab=expression(paste(delta^{13}, "C (‰, V-PDB)")))
legend("topleft", inset=c(0.05,0.0), legend=size.labels, col=colors, pch=19, cex=0.6, bty="n", x.intersp=1.8, y.intersp = 1.4)
arrows(mix.data$d13C-mix.data$d13C.se, mix.data$d15N, mix.data$d13C+mix.data$d13C.se, mix.data$d15N, length=0.0, angle=90, code=3, lwd=0.7, col=colors)
arrows(mix.data$d13C, mix.data$d15N-mix.data$d15N.se, mix.data$d13C, mix.data$d15N+mix.data$d15N.se, length=0.0, angle=90, code=3, lwd=0.7, col=colors)
points(d15N~d13C, data=mix.data, pch=19, cex=0.8, col=colors)
Figure Isotope values for size-fractioned seawater particles and plankters pooled across seasons and reef sites
dev.copy(pdf, "figures/environmental/iso.sources.KBay.pdf", encod="MacRoman", height=4, width=4)
dev.off()Biology
################################
### Biological responses
################################
#data : this is the master file
# add in light at depth column from df.light dataframe
data$Light<-df.light$Light
##### produce a categorical depth bin ####
depth<-data$newDepth
data$depth.bin<-factor(ifelse(depth<2, "<2m", ifelse(depth >2 & depth <4, "2-4m", ifelse(depth >4 & depth <6, "4-6m", ">6m"))), levels=c("<2m", "2-4m", "4-6m", ">6m"))
aggregate(Sample.ID~depth.bin+Season+Location, data, length)
data$depth.bin.small<-factor(ifelse(depth<4, "<4m", ">4m"), levels= c("<4m", ">4m"))
################################################
# calculate, normalized dependent variables
################################################
str(data)
data$cells.ml<-as.numeric(data$cells.ml)
# helpful shorthand
SA<-data$surface.area # surface area in cm2
blastate<-data$total.blastate.ml # tissue slurry blastate in ml
# AFDW.mg. == convert AFDW g to mg, mutiply by blastate volume, divide by cm2
data$biomass<- (data$mg.biomass.ml*blastate)/SA
# Symbiodinium.cells. == cell.ml * blastate / SA
data$zoox<- (data$cells.ml*blastate)/SA
# total chlorophyll == ug.chl.a.ml * blastate + ug.chl.c2.ml * blastate / SA
data$chltot<-(data$ug.chl.a.ml)+(data$ug.chl.c2.ml)*blastate/SA
# pg.chlorophyll.a..cell + pg.chlorophyll.c2..cell == ug.chltot.ml * 10^6 / cells.ml
data$chlcell<- (data$ug.chl.a.ml*10^6+data$ug.chl.c2.ml*10^6)/data$cells.mlqPCR
####################
# qPCR
#########
# qPCR
library(devtools)
# Use steponeR to import data and calculate proporation of C and D symbionts
source_url("https://raw.githubusercontent.com/jrcunning/steponeR/master/steponeR.R")
Mcap.plates <- list.files(path="data/qPCR", pattern = "txt$", full.names = T); Mcap.plates
Mcap <- steponeR(files=Mcap.plates, delim="\t",
target.ratios=c("C.D"),
fluor.norm=list(C=2.26827, D=0),
copy.number=list(C=33, D=3),
ploidy=list(C=1, D=1),
extract=list(C=0.813, D=0.813))
Mcap <- Mcap$result
head(Mcap)
# remove +/-control
Mcap <- Mcap[grep("+C52", Mcap$Sample.Name, fixed=T, invert = T), ]
Mcap <- Mcap[grep("H2O", Mcap$Sample.Name, fixed=T, invert = T), ]
# to remove any early-amplification CT noise
Mcap$C.CT.mean[which(Mcap$C.CT.mean < 15)] <- 0
#Remove failed samples, i.e., those where either C or D were NOT found in both reps
Mcap$fail <- ifelse(Mcap$C.reps < 2 & Mcap$D.reps < 2, TRUE, FALSE)
fails <- Mcap[Mcap$fail==TRUE, ]
Mcap <- Mcap[which(Mcap$fail==FALSE),]
# replace CT means with 'NA' as zero
Mcap$C.CT.mean[is.na(Mcap$C.CT.mean)] <-0
Mcap$D.CT.mean[is.na(Mcap$D.CT.mean)] <-0
Mcap$C.D[is.na(Mcap$C.D)] <- 1 # sets all infinity (= 100% C) to 1.0
# caluclate proportion C and proprtion D where C and D are both present
Mcap$propC<- Mcap$C.D / (Mcap$C.D + 1)
Mcap$propD<- 1 / (Mcap$C.D + 1)
# where C and D are not cooccuring...
# if C.D = 1 = 100% C, make 'PropC' = 1 and 'PropD' = 0
# if C.D = 0 = 100% D, make 'PropD' = 1 and 'PropC' = 0
Mcap$propC[which(Mcap$C.D==1)] <- 1
Mcap$propD[which(Mcap$propC==1)] <- 0
Mcap$propD[which(Mcap$C.D==0)] <- 1
# calculate FOUR COMMUNITY categories: C, C>D, D>C, D
Mcap$Mix <- factor(ifelse(Mcap$propC > Mcap$propD, ifelse(Mcap$propD!= 0, "CD", "C"), ifelse(Mcap$propD > Mcap$propC, ifelse(Mcap$propC!=0, "DC", "D"), NA)), levels=c("C", "CD", "DC", "D"))
# Identify SINGLE dominant symbiont clade: C or D
Mcap$Dom <- factor(substr(as.character(Mcap$Mix), 1, 1))
# Set zeros to NA to facilitate log transformation
Mcap$propC[which(Mcap$propC==0)] <- NA
Mcap$propD[which(Mcap$propD==0)] <- NA
######## look for duplicates in dataset by year and type of event (bleach/recover)
Mcap[duplicated(Mcap$Sample.Name), ] ## duplicates
# remove duplicated
Mcap<-Mcap[!(Mcap$Sample.Name=="HIMB_15" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),]
Mcap<-Mcap[!(Mcap$Sample.Name=="R10_05" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),]
Mcap<-Mcap[!(Mcap$Sample.Name=="R42_06" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),]
Mcap<-Mcap[!(Mcap$Sample.Name=="R46_01" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),]
Mcap<-Mcap[!(Mcap$Sample.Name=="R46_02" & Mcap$File.Name=="Wall_PanKbay_plate2.txt"),]
Mcap<-Mcap[!(Mcap$Sample.Name=="R46_03" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),]
Mcap<-Mcap[!(Mcap$Sample.Name=="HIMB_13" & Mcap$File.Name=="Wall_PanKbay_plate2.txt"),]
Mcap<-Mcap[!(Mcap$Sample.Name=="HIMB_14" & Mcap$File.Name=="Wall_PanKbay_plate2.txt"),]
Mcap<-Mcap[!(Mcap$Sample.Name=="R10_15" & Mcap$File.Name=="Wall_PanKbay_plate3.txt"),]
Mcap<-Mcap[!(Mcap$Sample.Name=="R10_15" & Mcap$File.Name=="Wall_PanKbay_plate1.txt"),]
Mcap<-Mcap[!(Mcap$Sample.Name=="R42_11" & Mcap$File.Name=="Wall_PanKbay_plate2.txt"),]
# parse Sample ID and Site from "Site.Name"
Mcap<-cbind(Mcap, colsplit(Mcap$Sample.Name, pattern= "_", c("Reef", "Sample.ID")))
Mcap$Season<-as.factor("winter")
Mcap$Reef<-as.factor(Mcap$Reef)
Mcap$Reef<-revalue(Mcap$Reef, c("R10"="F8-R10", "R46"="F1-R46")) # rename factor levels
# make new factors for bay region and reef type
Mcap$Bay.region <- ifelse(Mcap$Reef=="R42" | Mcap$Reef=="F1-R46", "northern", "southern")
Mcap$Reef.type <- ifelse(Mcap$Reef=="R42" | Mcap$Reef=="HIMB", "patch", "fringe")
Mcap$Location <- ifelse(Mcap$Reef=="F1-R46", "NW", ifelse(Mcap$Reef=="R42", "NE",
ifelse(Mcap$Reef=="F8-R10", "SW", "SE")))
### reorder columns and finish
Mcap<-Mcap[, c(17,15,20,19,18,16,1:14)] # reordered to match masterdata, and finish
### structure winter and summer qPCR dataframes to have same columns and combine dataframe
qPCR.winter<-Mcap[ , (names(Mcap) %in%
c("Season", "Reef", "Location", "Reef.type", "Bay.region", "Sample.ID", "propC", "propD", "Mix", "Dom"))]
qPCR.summer<-qPCR.Innis[ , (names(qPCR.Innis) %in%
c("Season", "Reef", "Location", "Reef.type", "Bay.region", "Sample.ID", "propC", "propD", "Mix", "Dom"))]
# merge qPCR files
qPCR.all<-rbind(qPCR.winter, qPCR.summer)
# add to master data
data.all<-merge(data, qPCR.all, by=c("Season", "Reef", "Location", "Reef.type", "Bay.region", "Sample.ID"), all.x=T)
###### remove columns no longer needed, update "Depth" to be tide-corrected depth )= newDepth
data.trim<-data.all[ , !(names(data.all) %in% c("total.blastate.ml", "Date", "Time.of.collection", "Depth..m", "Time.r", "surface.area.cm2", "cells.ml", "ug.chl.a.ml", "ug.chl.c2.ml", "mg.biomass.ml", "host..mass.mg", "host..ugN", "host..ugC", "symb..mass.mg", "symb..ugN", "symb..ugC", "stationId", "datum", "TimeUTC", "TideHT", "Time"))]
data.trim$symb..C.N[data.trim$symb..C.N>=12.520270]=NA # set this outlier to NAqPCR figure The distribution of C and D symbiont types is modeled using a logistic regression with dominant symbiont community (C vs. D) as a response and Depth (m), Season, and reef Location as a predictor. A generalized linear model with a binomial distribution and logit link function is used. Model selection (by model AIC values) between a fully crossed model (crossed) and additive model (add) were compared. The additive mode best fit the data; main effects were tested in anova with a Chi-square test.
## qPCR figures
#####
model.data<-data.trim[,c(17, 1:5, 18:25, 7:16, 26:29)]
#########################
#########################
#########################
# Plots with DEPTH as the predictor
# qPCR and symbionts
#########################
###### Plot Dominant Symbiont and Depth (both seasons)
Symb<-model.data
Symb$Reef <- factor(Symb$Reef, levels = c("F1-R46", "R42", "F8-R10", "HIMB")) # set levels
Symb$Location <- factor(Symb$Location, levels = c("NW", "NE", "SW", "SE")) # set levels
Symb$Dominant <- ifelse(Symb$Dom=="D", 1, 0)
add <-glm(Dominant~newDepth+Season+Location, family = "binomial", data = Symb)
crossed <-glm(Dominant~newDepth*Season+Location, family = "binomial", data = Symb)
AIC(add, crossed) # additive model best## df AIC
## add 6 129.7125
## crossed 7 125.5697anova(crossed, test = "Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Dominant
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 119 143.06
## newDepth 1 21.6786 118 121.38 3.224e-06 ***
## Season 1 0.3020 117 121.08 0.5826
## Location 3 3.3647 114 117.71 0.3387
## newDepth:Season 1 6.1428 113 111.57 0.0132 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1results.all<-crossed
#summary(results.all)
plot(allEffects(results.all), par.strip.text=list(cex=0.6))
results.all.mod <-glm(Dominant~newDepth, family = "binomial", data = Symb)
fitted.all <- predict(results.all.mod, newdata = list(newDepth=seq(0,10,0.1)), type = "response")
###### summer only symbionts
sum.dat<-Symb[(Symb$Season=="summer"),]
results.sum=glm(Dominant~newDepth, family = "binomial", data = sum.dat)
anova(results.sum, test = "Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Dominant
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 59 67.48
## newDepth 1 21.71 58 45.77 3.171e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#summary(results.sum)
fitted.sum <- predict(results.sum, newdata = list(newDepth=seq(0,10,0.1)), type = "response")
# plot(fitted.sum, ylab="proportion D", ylim=c(0,1))
###### winter only symbionts
wint.dat<-Symb[(Symb$Season=="winter"),]
results.win=glm(Dominant~newDepth, family = "binomial", data = wint.dat)
anova(results.win, test = "Chisq")## Analysis of Deviance Table
##
## Model: binomial, link: logit
##
## Response: Dominant
##
## Terms added sequentially (first to last)
##
##
## Df Deviance Resid. Df Resid. Dev Pr(>Chi)
## NULL 59 74.920
## newDepth 1 4.6265 58 70.293 0.03148 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#summary(results.win)
fitted.win <- predict(results.win, newdata = list(newDepth=seq(0,10,0.1)), type = "response")
# plot(fitted.win, ylab="proportion D", ylim=c(0,1))Logistic regression of the probability of Montipora capitata have a clade D-dominant symbiont community (probability of 1.0) versus one dominated by clade C symbiont (probability of 0). Both seasons show similar patterns, although the probability of corals being clade D-dominant increases in winter months as a function of Depth relative to summer months.
######
######
## Figure of Dominant symbiont clades across seasons
## **Note** where points equal 0.0 is 100% D, where they equal 0 is 100% C
par(mar=c(5,4,3,2))
plot(sum.dat$propD~sum.dat$newDepth, xlab="Depth (m)", ylab = "Proportion of Clade D Symbiont", pch=19, col="salmon", xlim=c(0,10), ylim=c(0.0, 1.0), cex=0.5, cex.lab=0.8)
par(new=T)
plot(wint.dat$propD~wint.dat$newDepth, pch=19, col="lightskyblue", xlim=c(0,10), ylim=c(0.0, 1.0), xaxt="n", xlab="", ylab="", cex=0.5, cex.lab=0.8)
lines(fitted.all ~ seq(0,10,0.1), col="gray30", lwd=1, lty=2)
lines(fitted.sum ~ seq(0,10,0.1), col="coral", lwd=1)
lines(fitted.win ~ seq(0,10,0.1), col="lightskyblue", lwd=1)
legend("topright", pch=c(19,19, NA), lty=c(1,1,2), col=c("coral", "lightskyblue", "gray30"), legend=c("Summer", "Winter", "Combined"), lwd=1, bty="n", x.intersp=1, pt.cex=0.8, y.intersp=2, cex=0.5, inset=c(-0.005, -0.01), seg.len=4)
Figure Probability of Montipora capitata hosting clade D symbionts across a depth gradient in summer and winter seasons. Lines represent model fit to data by logistic regression for summer alone, winter alone, and combined (summer + winter) data
dev.copy(pdf, "figures/symbionts/Symbionts_by_Season_depth.pdf", encod="MacRoman", height=4, width=4)
dev.off()Seen another way–Histograms. By separating these data into 2 histograms the probability of clade D- (1.0) and clade C-dominance (0.0) can be visualized. Here it is clear that there is a greater probability of finding corals with clade D dominant communites at shallow depths, whereas in winter there is a slightly greater probability of finding clade D in corals with increasing depth.
#######
#######
par(mfrow=c(1,2))
Dom <- subset(Symb, !is.na(newDepth) & !is.na(Dominant))
Dom.sum<-Dom[(Dom$Season=="summer"),]
Dom.win<-Dom[(Dom$Season=="winter"),]
logi.hist.plot(Dom.sum$newDepth, Dom.sum$Dominant, boxp = FALSE, type = "hist", col="coral", xlabel = "Depth (m)", ylabel = "", ylabel2 = "", main="summer")
mtext(side = 2, text = "Probability of Clade D", line = 3, cex = 1)
logi.hist.plot(Dom.win$newDepth, Dom.win$Dominant, boxp = FALSE, type = "hist", col="lightskyblue", xlabel = "Depth (m)", ylabel = "", ylabel2 = "", main="winter")
mtext(side = 4, text = "Frequency", line = 0.5, cex=1)
mtext(side = 4, text = "C D", line = 0.5, cex = 0.8)
Figure Histograms of probability of Montipora capitata hosting clade D symbionts across a depth gradient among seasons. Lines represent model fit to data by logistic regression
dev.copy(pdf, "figures/symbionts/Symb_Season_Light_logistic.pdf", encod="MacRoman", height=5, width=8)
dev.off()#
#########################
#########################
# this code runs the same as above, except using colony LIGHT-AT-DEPTH as a predictor instead of depth
# Depth as a predictor
# qPCR and symbionts
#########################
###### Plot Dominant Symbiont and Depth (both seasons)
Symb<-model.data
Symb$Dominant <- ifelse(Symb$Dom=="D", 1, 0)
add<-glm(Dominant~Light+Season+Location, family = "binomial", data = Symb)
crossed<-glm(Dominant~Light*Season*Location, family = "binomial", data = Symb)
#AIC(add, crossed)
anova(add, test = "Chisq")
#summary(add)
results.all<-add
results.all.mod<-glm(Dominant~Light, family = "binomial", data = Symb)
fitted.all <- predict(results.all.mod, newdata = list(Light=seq(0,25,0.1)), type = "response")
###### summer only symbionts
sum.dat<-Symb[(Symb$Season=="summer"),]
results.sum<-glm(Dominant~Light, family = "binomial", data = sum.dat)
#anova(results.sum, test = "Chisq")
#summary(results.sum)
fitted.sum <- predict(results.sum, newdata = list(Light=seq(0,25,0.1)), type = "response")
# plot(fitted.sum, ylab="proportion D", ylim=c(0,1))
###### winter only symbionts
wint.dat<-Symb[(Symb$Season=="winter"),]
results.win<-glm(Dominant~Light, family = "binomial", data = wint.dat)
#anova(results.win, test = "Chisq")
#summary(results.win)
fitted.win <- predict(results.win, newdata = list(Light=seq(0,25,0.1)), type = "response")
# plot(fitted.win, ylab="proportion D", ylim=c(0,1))
######
######
## Figure of Dominant symbiont clades across seasons
## **Note** where points equal 0.0 is 100% D, where they equal 0 is 100% C
par(mar=c(5,4,3,2))
plot(sum.dat$propD~sum.dat$Light, xlab="Light (DLI)", ylab = "Proportion of Clade D Symbiont", pch=19, col="salmon", xlim=c(0,25), ylim=c(0.0, 1.0), cex=0.5, cex.lab=0.8)
par(new=T)
plot(wint.dat$propD~wint.dat$newDepth, pch=19, col="lightskyblue", xlim=c(0,25), ylim=c(0.0, 1.0), xaxt="n", xlab="", ylab="", cex=0.5)
lines(fitted.all ~ seq(0,25,0.1), col="gray30", lwd=1, lty=2)
lines(fitted.sum ~ seq(0,25,0.1), col="coral", lwd=1)
lines(fitted.win ~ seq(0,25,0.1), col="lightskyblue", lwd=1)
legend("bottomright", pch=c(19,19, NA), lty=c(1,1,2), col=c("coral", "lightskyblue", "gray30"), legend=c("Summer", "Winter", "Combined"), lwd=1, bty="n", x.intersp=1, pt.cex=0.7, y.intersp=2, cex=0.5, inset=c(0.1, 0.15), seg.len=4)
dev.copy(pdf, "figures/symbionts/Symbionts_by_Season_light.pdf", encod="MacRoman", height=4, width=4)
dev.off()
#######
#######
par(mfrow=c(1,2))
Dom <- subset(Symb, !is.na(Light) & !is.na(Dominant))
Dom.sum<-Dom[(Dom$Season=="summer"),]
Dom.win<-Dom[(Dom$Season=="winter"),]
logi.hist.plot(Dom.sum$Light, Dom.sum$Dominant, boxp = FALSE, type = "hist", col="coral", xlabel = "Light-at-depth (DLI)", ylabel = "", ylabel2 = "", main="summer")
mtext(side = 2, text = "Probability of Clade D", line = 3, cex = 1)
logi.hist.plot(Dom.win$Light, Dom.win$Dominant, boxp = FALSE, type = "hist", col="lightskyblue", xlabel = "Light-at-depth (DLI)", ylabel = "", ylabel2 = "", main="winter")
mtext(side = 4, text = "Frequency", line = 0.5, cex=1)
mtext(side = 4, text = "C D", line = 0.5, cex = 0.8)
dev.copy(pdf, "figures/symbionts/Symb_Season_Light_logistic.pdf", encod="MacRoman", height=5, width=8)
dev.off()Model assumptions
Physiology models
Total Biomass
########## ##########
######### biomass ----
Y<-model.data$biomass
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use additive model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 6 846.22 862.95 -417.11 834.22
## object 10 848.36 876.23 -414.18 828.36 5.8663 4 0.2094summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 828.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.12288 -0.72960 -0.07959 0.55985 3.08964
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 5.82 2.412
## Residual 60.57 7.783
## Number of obs: 120, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 25.62373 2.24200 16.01000 11.429 4.14e-09 ***
## Seasonwinter 0.08287 1.61545 115.55000 0.051 0.959
## Light -0.04636 0.17286 108.20000 -0.268 0.789
## DomD 2.87278 1.75150 115.80000 1.640 0.104
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light
## Seasonwintr -0.593
## Light -0.692 0.471
## DomD 0.131 -0.258 -0.421rand(add)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## Location 4.22 1 0.04 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(anova(add, type=2), digits=5)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 0.159 0.159 1 115.55 0.00263 0.9592
## Light 4.358 4.358 1 108.20 0.07194 0.7890
## Dom 162.958 162.958 1 115.80 2.69019 0.1037var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.biomass=var.table
var.biomass$response<-"biomass"
plot(allEffects(add), ylab="biomass", par.strip.text=list(cex=0.7))
Symbiont density
######### symbionts--
Y<-model.data$zoox
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+Season:Light + (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use additive model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + Season:Light + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 7 16.414 35.926 -1.20681 2.4136
## object 10 21.951 49.826 -0.97577 1.9515 0.4621 3 0.9271summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + Season:Light + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 30.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.08854 -0.62190 -0.06602 0.57425 2.14357
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.02478 0.1574
## Residual 0.05715 0.2391
## Number of obs: 120, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 14.309970 0.101529 6.920000 140.945 3.19e-13 ***
## Seasonwinter -0.031559 0.086562 112.270000 -0.365 0.71611
## Light 0.005428 0.005991 113.810000 0.906 0.36688
## DomD 0.401155 0.054175 113.140000 7.405 2.52e-11 ***
## Seasonwinter:Light 0.029607 0.010208 113.290000 2.900 0.00448 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light DomD
## Seasonwintr -0.439
## Light -0.540 0.593
## DomD 0.099 -0.153 -0.393
## Ssnwntr:Lgh 0.247 -0.817 -0.417 -0.001rand(add)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## Location 22.8 1 2e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(anova(add, type=2), digits=5)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 0.00760 0.00760 1 112.27 0.133 0.716111
## Light 0.39237 0.39237 1 114.92 6.866 0.009971 **
## Dom 3.13331 3.13331 1 113.14 54.830 2.521e-11 ***
## Season:Light 0.48067 0.48067 1 113.29 8.411 0.004480 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.symb=var.table
var.symb$response<-"symb"
posthoc<-emmeans(add, ~Dom)
CLD(posthoc, Letters=letters)## Dom emmean SE df lower.CL upper.CL .group
## C 14.45561 0.08380064 3.26 14.20043 14.71080 a
## D 14.85677 0.09108456 4.53 14.61520 15.09834 b
##
## Results are averaged over the levels of: Season
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(add), ylab="symbiont density", par.strip.text=list(cex=0.7))
Chlorophyll a
######### chltotal--
Y<-model.data$chltot
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+Season:Dom + (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use additive model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + Season:Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 7 424.92 444.44 -205.46 410.92
## object 10 425.73 453.61 -202.87 405.73 5.1938 3 0.1581summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + Season:Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 418.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.04486 -0.63880 -0.07154 0.47425 2.98104
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.2169 0.4657
## Residual 1.7844 1.3358
## Number of obs: 120, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 5.84785 0.40030 13.08000 14.609 1.77e-09 ***
## Seasonwinter 1.70590 0.30499 113.71000 5.593 1.56e-07 ***
## Light -0.09232 0.03079 109.93000 -2.998 0.00336 **
## DomD -0.49207 0.44865 114.17000 -1.097 0.27504
## Seasonwinter:DomD -1.32504 0.56779 113.45000 -2.334 0.02137 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light DomD
## Seasonwintr -0.528
## Light -0.644 0.315
## DomD 0.074 0.148 -0.459
## Ssnwntr:DmD 0.017 -0.414 0.247 -0.741rand(add)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## Location 4.87 1 0.03 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(anova(add, type=2), digits=5)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 45.581 45.581 1 114.27 25.5446 1.658e-06 ***
## Light 16.042 16.042 1 109.93 8.9901 0.003357 **
## Dom 31.622 31.622 1 114.55 17.7213 5.111e-05 ***
## Season:Dom 9.718 9.718 1 113.45 5.4461 0.021372 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.chl=var.table
var.chl$response<-"chl"
posthoc<-emmeans(add, ~Season)
CLD(posthoc, Letters=letters)## Season emmean SE df lower.CL upper.CL .group
## summer 4.865228 0.3288109 6.70 4.080513 5.649944 a
## winter 5.908608 0.3026093 4.98 5.129885 6.687331 b
##
## Results are averaged over the levels of: Dom
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05posthoc<-emmeans(add, ~Dom|Season)
CLD(posthoc, Letters=letters)## Season = summer:
## Dom emmean SE df lower.CL upper.CL .group
## D 4.619195 0.4738171 24.45 3.642243 5.596147 a
## C 5.111262 0.3066745 5.30 4.336006 5.886518 a
##
## Season = winter:
## Dom emmean SE df lower.CL upper.CL .group
## D 5.000055 0.3888373 12.97 4.159808 5.840303 a
## C 6.817161 0.3257516 6.50 6.034599 7.599723 b
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(add), ylab="total chlor", par.strip.text=list(cex=0.7))
Chlorophyll per symbiont cell
######### chlcell --
Y<-model.data$chlcell
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use additive model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 6 4.3051 21.030 3.8475 -7.6949
## object 10 6.3305 34.205 6.8347 -13.6695 5.9745 4 0.2011summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 13.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.73631 -0.52024 -0.08737 0.53912 2.46089
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.02493 0.1579
## Residual 0.05191 0.2278
## Number of obs: 120, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.874743 0.096951 5.720000 19.337 1.98e-06 ***
## Seasonwinter 0.079964 0.047578 113.860000 1.681 0.0956 .
## Light -0.035217 0.005194 115.720000 -6.781 5.33e-10 ***
## DomD -0.513165 0.051649 114.030000 -9.936 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light
## Seasonwintr -0.411
## Light -0.480 0.481
## DomD 0.100 -0.267 -0.433rand(add)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## Location 23.7 1 1e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(anova(add, type=2), digits=5)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 0.1466 0.1466 1 113.86 2.825 0.09556 .
## Light 2.3869 2.3869 1 115.72 45.977 5.331e-10 ***
## Dom 5.1248 5.1248 1 114.03 98.716 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.chlcell=var.table
var.chlcell$response<-"chlcell"
plot(allEffects(add), ylab="chla cell", par.strip.text=list(cex=0.7))
#ranef(add)
#rand(add)
#fixef(add)
#sjp.lmer(add, y.offset = .4)Isotope models
host d13C
########## host..d13C --
Y<-model.data$host..d13C
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ Season:Dom+ Season:Light+(1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use additive model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + Season:Dom + Season:Light + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 8 323.78 346.08 -153.89 307.78
## object 10 327.70 355.58 -153.85 307.70 0.0775 2 0.962summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula:
## Y ~ Season + Light + Dom + Season:Dom + Season:Light + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 323.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.77180 -0.68022 -0.08845 0.75170 2.22371
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.2658 0.5156
## Residual 0.7386 0.8594
## Number of obs: 120, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -16.41116 0.34694 7.99000 -47.303 4.47e-11 ***
## Seasonwinter 0.01219 0.31119 111.36000 0.039 0.9688
## Light 0.13233 0.02315 112.82000 5.717 9.01e-08 ***
## DomD -1.62857 0.29962 111.59000 -5.436 3.25e-07 ***
## Seasonwinter:DomD 0.84855 0.38854 111.43000 2.184 0.0311 *
## Seasonwinter:Light -0.01870 0.03892 112.30000 -0.480 0.6319
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light DomD Ssn:DD
## Seasonwintr -0.463
## Light -0.558 0.557
## DomD 0.130 -0.110 -0.517
## Ssnwntr:DmD -0.083 0.014 0.368 -0.760
## Ssnwntr:Lgh 0.274 -0.775 -0.491 0.254 -0.335rand(add)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## Location 20.1 1 7e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(anova(add), digits=5)## Analysis of Variance Table of type III with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 1.0324 1.0324 1 111.22 1.398 0.23961
## Light 23.6384 23.6384 1 113.94 32.005 1.160e-07 ***
## Dom 27.4555 27.4555 1 112.24 37.174 1.566e-08 ***
## Season:Dom 3.5227 3.5227 1 111.43 4.770 0.03106 *
## Season:Light 0.1704 0.1704 1 112.30 0.231 0.63192
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.d13C.H=var.table
var.d13C.H$response<-"d13C.H"
posthoc<-emmeans(add, ~Dom)
CLD(posthoc, Letters=letters)## Dom emmean SE df lower.CL upper.CL .group
## D -16.62814 0.3090539 5.04 -17.42049 -15.83579 a
## C -15.42384 0.2779424 3.31 -16.26338 -14.58431 b
##
## Results are averaged over the levels of: Season
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05posthoc<-emmeans(add, ~Dom| Season)
CLD(posthoc, Letters=letters)## Season = summer:
## Dom emmean SE df lower.CL upper.CL .group
## D -16.98392 0.3780843 10.98 -17.81631 -16.15154 a
## C -15.35535 0.2880146 3.83 -16.16948 -14.54122 b
##
## Season = winter:
## Dom emmean SE df lower.CL upper.CL .group
## D -16.27236 0.3270155 6.30 -17.06346 -15.48125 a
## C -15.49233 0.3058112 4.79 -16.28894 -14.69573 b
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(add), ylab="d13Chost", par.strip.text=list(cex=0.7))
symbiont d13C
########## symb..d13C --
Y<-model.data$symb..d13C
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ Season:Dom +(1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use additive model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + Season:Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 7 334.15 353.66 -160.08 320.15
## object 10 339.80 367.68 -159.90 319.80 0.3464 3 0.9511summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + Season:Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 330.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9724 -0.6051 -0.0345 0.7590 2.8420
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.3887 0.6234
## Residual 0.8043 0.8969
## Number of obs: 120, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -16.59013 0.38242 5.80000 -43.382 1.64e-08 ***
## Seasonwinter -0.29714 0.20544 112.63000 -1.446 0.1509
## Light 0.14086 0.02111 114.74000 6.673 9.30e-10 ***
## DomD -1.53545 0.30257 112.84000 -5.075 1.54e-06 ***
## Seasonwinter:DomD 1.26115 0.38222 112.51000 3.300 0.0013 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light DomD
## Seasonwintr -0.377
## Light -0.462 0.322
## DomD 0.061 0.141 -0.466
## Ssnwntr:DmD 0.009 -0.411 0.249 -0.741rand(add)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## Location 26.9 1 2e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(anova(add, type=2), digits=5)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 0.002 0.002 1 112.90 0.002 0.9615036
## Light 35.816 35.816 1 114.74 44.529 9.298e-10 ***
## Dom 12.375 12.375 1 113.05 15.386 0.0001508 ***
## Season:Dom 8.757 8.757 1 112.51 10.887 0.0012975 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.d13C.S=var.table
var.d13C.S$response<-"d13C.S"
posthoc<-emmeans(add, ~Dom)
CLD(posthoc, Letters=letters)## Dom emmean SE df lower.CL upper.CL .group
## D -16.51965 0.3551755 4.34 -17.47613 -15.56317 a
## C -15.61478 0.3274383 3.15 -16.62935 -14.60020 b
##
## Results are averaged over the levels of: Season
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05posthoc<-emmeans(add, ~Dom| Season)
CLD(posthoc, Letters=letters)## Season = summer:
## Dom emmean SE df lower.CL upper.CL .group
## D -17.00166 0.4167784 8.07 -17.96137 -16.04195 a
## C -15.46621 0.3392889 3.63 -16.44701 -14.48541 b
##
## Season = winter:
## Dom emmean SE df lower.CL upper.CL .group
## D -16.03764 0.3752434 5.40 -16.98092 -15.09436 a
## C -15.76335 0.3471338 3.96 -16.73061 -14.79608 a
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(add), ylab="d13Csymb", par.strip.text=list(cex=0.7))
host-symbiont d13C
######### d13C..host.sym --
Y<-model.data$d13C..host.sym
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ Season:Light +Season:Dom +(1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use add model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + Season:Light + Season:Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 8 116.04 138.34 -50.022 100.044
## object 10 119.77 147.64 -49.884 99.768 0.2757 2 0.8712summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula:
## Y ~ Season + Light + Dom + Season:Light + Season:Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 126.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.9375 -0.3451 0.0104 0.4548 3.6783
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.02824 0.1680
## Residual 0.13287 0.3645
## Number of obs: 120, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.126970 0.129306 10.970000 0.982 0.34730
## Seasonwinter 0.409353 0.131934 111.360000 3.103 0.00243 **
## Light -0.002278 0.009790 113.440000 -0.233 0.81639
## DomD -0.136399 0.126986 111.720000 -1.074 0.28508
## Seasonwinter:Light -0.034256 0.016478 112.780000 -2.079 0.03990 *
## Seasonwinter:DomD -0.347111 0.164712 111.490000 -2.107 0.03733 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light DomD Ssnw:L
## Seasonwintr -0.528
## Light -0.633 0.560
## DomD 0.147 -0.111 -0.516
## Ssnwntr:Lgh 0.315 -0.775 -0.497 0.257
## Ssnwntr:DmD -0.094 0.015 0.370 -0.761 -0.336rand(add)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## Location 8.21 1 0.004 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(anova(add, type=2), digits=4)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 1.3196 1.3196 1 111.4 9.931 0.00209 **
## Light 0.3604 0.3604 1 112.9 2.712 0.10235
## Dom 2.2911 2.2911 1 112.8 17.243 6.42e-05 ***
## Season:Light 0.5742 0.5742 1 112.8 4.322 0.03990 *
## Season:Dom 0.5901 0.5901 1 111.5 4.441 0.03733 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.d13C.HS=var.table
var.d13C.HS$response<-"d13C.HS"
posthoc<-emmeans(add, ~Dom)
CLD(posthoc, Letters=letters)## Dom emmean SE df lower.CL upper.CL .group
## D -0.1331473 0.11086251 6.46 -0.3998082 0.1335136 a
## C 0.1768072 0.09489086 3.49 -0.1025094 0.4561238 b
##
## Results are averaged over the levels of: Season
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05posthoc<-emmeans(add, ~Dom|Season)
CLD(posthoc, Letters=letters)## Season = summer:
## Dom emmean SE df lower.CL upper.CL .group
## D -0.02760769 0.1443214 17.32 -0.33166817 0.27645279 a
## C 0.10879151 0.1001424 4.36 -0.16041776 0.37800078 a
##
## Season = winter:
## Dom emmean SE df lower.CL upper.CL .group
## D -0.23868691 0.1197827 8.71 -0.51101369 0.03363987 a
## C 0.24482291 0.1092528 5.96 -0.02294213 0.51258796 b
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(add), ylab="d13C-hs", par.strip.text=list(cex=0.7))
skeleton d13C
####### d13C..skel --
Y<-model.data$d13C..skel
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use add model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 6 313.69 330.42 -150.85 301.69
## object 10 319.03 346.91 -149.52 299.03 2.6592 4 0.6164summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 313
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.96117 -0.81100 -0.02084 0.65228 2.40914
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.1400 0.3742
## Residual 0.7022 0.8380
## Number of obs: 120, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -2.855001 0.277939 9.890000 -10.272 1.36e-06 ***
## Seasonwinter 0.460523 0.174545 114.750000 2.638 0.00949 **
## Light 0.008887 0.018898 115.180000 0.470 0.63906
## DomD 0.011169 0.189383 115.010000 0.059 0.95308
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light
## Seasonwintr -0.522
## Light -0.610 0.477
## DomD 0.122 -0.263 -0.428rand(add)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## Location 10.5 1 0.001 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(anova(add, type=2), digits=4)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 4.888 4.888 1 114.8 6.961 0.00949 **
## Light 0.155 0.155 1 115.2 0.221 0.63906
## Dom 0.002 0.002 1 115.0 0.003 0.95308
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.d13C.sk=var.table
var.d13C.sk$response<-"d13C.sk"
posthoc<-emmeans(add, ~Season)
CLD(posthoc, Letters=letters)## Season emmean SE df lower.CL upper.CL .group
## summer -2.778508 0.2286013 4.69 -3.377847 -2.179170 a
## winter -2.317985 0.2197840 4.07 -2.924134 -1.711836 b
##
## Results are averaged over the levels of: Dom
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(add), ylab="d13C-skel", par.strip.text=list(cex=0.7))
host d15N
####### host..d15N --
Y<-model.data$host..d15N
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+(1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use add model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 6 86.080 102.81 -37.040 74.080
## object 10 88.549 116.42 -34.275 68.549 5.5304 4 0.2371summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 90.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.4078 -0.5566 0.0705 0.6116 2.4030
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.29199 0.5404
## Residual 0.09669 0.3110
## Number of obs: 120, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.867242 0.280992 3.430000 17.322 0.000186 ***
## Seasonwinter 0.069196 0.065050 113.160000 1.064 0.289709
## Light -0.014857 0.007142 113.660000 -2.080 0.039755 *
## DomD 0.044236 0.070641 113.200000 0.626 0.532435
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light
## Seasonwintr -0.195
## Light -0.228 0.484
## DomD 0.049 -0.270 -0.436rand(add)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## Location 121 1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(anova(add, type=2), digits=4)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 0.1094 0.1094 1 113.2 1.132 0.2897
## Light 0.4184 0.4184 1 113.7 4.327 0.0398 *
## Dom 0.0379 0.0379 1 113.2 0.392 0.5324
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.d15N.H=var.table
var.d15N.H$response<-"d15N.H"
plot(allEffects(add), ylab="d15Nhost", par.strip.text=list(cex=0.7))
symbiont d15N
######## symb..d15N --
Y<-model.data$symb..d15N
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ Season:Dom+ Light:Dom + (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use add model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + Season:Dom + Light:Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 8 103.85 126.15 -43.925 87.850
## object 10 107.53 135.40 -43.763 87.526 0.3242 2 0.8504summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula:
## Y ~ Season + Light + Dom + Season:Dom + Light:Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 112.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.10105 -0.57497 0.06009 0.58736 2.67836
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.4296 0.6555
## Residual 0.1094 0.3307
## Number of obs: 120, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 4.343041 0.339414 3.390000 12.796 0.000551 ***
## Seasonwinter -0.080336 0.076829 111.090000 -1.046 0.297994
## Light -0.004405 0.008747 111.420000 -0.504 0.615491
## DomD -0.183898 0.247719 111.040000 -0.742 0.459434
## Seasonwinter:DomD 0.435748 0.179700 111.070000 2.425 0.016926 *
## Light:DomD 0.013734 0.016369 111.020000 0.839 0.403249
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light DomD Ssn:DD
## Seasonwintr -0.171
## Light -0.216 0.358
## DomD -0.073 0.203 0.204
## Ssnwntr:DmD 0.061 -0.416 -0.098 -0.815
## Light:DomD 0.095 -0.158 -0.442 -0.892 0.619rand(add)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## Location 139 1 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(anova(add, type=2), digits=6)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 0.000027 0.000027 1 111.116 0.00025 0.987387
## Light 0.002411 0.002411 1 111.518 0.02204 0.882246
## Dom 0.506556 0.506556 1 111.038 4.63087 0.033568 *
## Season:Dom 0.643192 0.643192 1 111.066 5.87997 0.016926 *
## Light:Dom 0.077007 0.077007 1 111.016 0.70399 0.403249
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.d15N.S=var.table
var.d15N.S$response<-"d15N.S"
posthoc<-emmeans(add, ~Dom|Season)
CLD(posthoc, Letters=letters)## Season = summer:
## Dom emmean SE df lower.CL upper.CL .group
## D 4.233574 0.3545527 4.03 3.251927 5.215221 a
## C 4.307890 0.3314409 3.08 3.268355 5.347426 a
##
## Season = winter:
## Dom emmean SE df lower.CL upper.CL .group
## C 4.227554 0.3327041 3.13 3.192580 5.262528 a
## D 4.588986 0.3368169 3.28 3.567657 5.610315 b
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(add), ylab="d15Nsymb", par.strip.text=list(cex=0.7))
host-symbiont d15N
####### d15N..host.sym --
Y<-model.data$d15N..host.sym
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ Season:Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use add model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + Season:Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 7 11.866 31.378 1.0670 -2.1341
## object 10 15.598 43.472 2.2013 -4.4025 2.2685 3 0.5186summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + Season:Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: 21.9
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.0812 -0.6452 0.1012 0.7207 2.0871
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.01312 0.1145
## Residual 0.05606 0.2368
## Number of obs: 120, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 0.538132 0.081629 8.560000 6.592 0.000126 ***
## Seasonwinter 0.162519 0.054163 113.070000 3.001 0.003316 **
## Light -0.013270 0.005526 114.610000 -2.401 0.017938 *
## DomD 0.078624 0.079731 113.430000 0.986 0.326175
## Seasonwinter:DomD -0.417705 0.100797 112.870000 -4.144 6.62e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light DomD
## Seasonwintr -0.463
## Light -0.566 0.319
## DomD 0.071 0.144 -0.463
## Ssnwntr:DmD 0.012 -0.412 0.248 -0.741rand(add)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## Location 10.6 1 0.001 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(anova(add, type=2), digits=4)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 0.1037 0.1037 1 113.5 1.849 0.17659
## Light 0.3233 0.3233 1 114.6 5.767 0.01794 *
## Dom 0.5375 0.5375 1 113.8 9.588 0.00247 **
## Season:Dom 0.9627 0.9627 1 112.9 17.173 6.62e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.d15N.HS<-var.table
var.d15N.HS$response<-"d15N.HS"
posthoc<-emmeans(add, ~Dom|Season)
CLD(posthoc, Letters=letters)## Season = summer:
## Dom emmean SE df lower.CL upper.CL .group
## C 0.4322534 0.06730944 4.27 0.24988899 0.6146178 a
## D 0.5108774 0.09284886 14.23 0.31203153 0.7097232 a
##
## Season = winter:
## Dom emmean SE df lower.CL upper.CL .group
## D 0.2556905 0.07951643 8.15 0.07290503 0.4384760 a
## C 0.5947720 0.07005404 4.94 0.41398875 0.7755552 b
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
## significance level used: alpha = 0.05plot(allEffects(add), ylab="d15N-hs", par.strip.text=list(cex=0.7))
host C:N
###### host..C.N --
Y<-model.data$host..C.N
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use add model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 6 -273.59 -256.87 142.79 -285.59
## object 10 -268.55 -240.68 144.28 -288.55 2.9622 4 0.5642summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: -254.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7879 -0.7136 -0.1026 0.6668 3.2204
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.0008847 0.02974
## Residual 0.0052870 0.07271
## Number of obs: 120, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1.944e+00 2.319e-02 1.050e+01 83.819 4.44e-16 ***
## Seasonwinter 3.999e-03 1.513e-02 1.149e+02 0.264 0.792
## Light 1.371e-03 1.635e-03 1.142e+02 0.839 0.403
## DomD 6.456e-03 1.642e-02 1.152e+02 0.393 0.695
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light
## Seasonwintr -0.541
## Light -0.632 0.476
## DomD 0.125 -0.262 -0.427rand(add)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## Location 6.87 1 0.009 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1print(anova(add, type=2), digits=4)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 0.000369 0.000369 1 114.9 0.0698 0.792
## Light 0.003718 0.003718 1 114.2 0.7032 0.403
## Dom 0.000818 0.000818 1 115.2 0.1546 0.695var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.CN.H<-var.table
var.CN.H$response<-"CN.Host"
plot(allEffects(add), ylab="C:N host", par.strip.text=list(cex=0.7))
symbiont C:N
####### symb..C.N --
Y<-model.data$symb..C.N
full<-lmer(Y~Season*Light*Dom+ (1|Location), data=model.data, na.action=na.exclude)
add<-lmer(Y~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
anova(full, add) #use add model## Data: model.data
## Models:
## ..1: Y ~ Season + Light + Dom + (1 | Location)
## object: Y ~ Season * Light * Dom + (1 | Location)
## Df AIC BIC logLik deviance Chisq Chi Df Pr(>Chisq)
## ..1 6 -212.29 -195.61 112.14 -224.29
## object 10 -207.16 -179.37 113.58 -227.16 2.8707 4 0.5797summary(add)## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: Y ~ Season + Light + Dom + (1 | Location)
## Data: model.data
##
## REML criterion at convergence: -193.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.8398 -0.6571 -0.1712 0.6883 2.9658
##
## Random effects:
## Groups Name Variance Std.Dev.
## Location (Intercept) 0.000000 0.00000
## Residual 0.009201 0.09592
## Number of obs: 119, groups: Location, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.051e+00 2.207e-02 1.150e+02 92.903 <2e-16 ***
## Seasonwinter -2.956e-02 1.957e-02 1.150e+02 -1.510 0.134
## Light 3.778e-04 1.953e-03 1.150e+02 0.193 0.847
## DomD 8.659e-04 2.110e-02 1.150e+02 0.041 0.967
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) Ssnwnt Light
## Seasonwintr -0.684
## Light -0.798 0.433
## DomD 0.097 -0.230 -0.379rand(add)## Analysis of Random effects Table:
## Chi.sq Chi.DF p.value
## Location -2.84e-14 1 1print(anova(add, type=2), digits=5)## Analysis of Variance Table of type II with Satterthwaite
## approximation for degrees of freedom
## Sum Sq Mean Sq NumDF DenDF F.value Pr(>F)
## Season 0.0209902 0.0209902 1 115 2.28128 0.1337
## Light 0.0003443 0.0003443 1 115 0.03742 0.8470
## Dom 0.0000155 0.0000155 1 115 0.00168 0.9673var.table<-as.data.frame(VarCorr(add),comp="Variance")
var.table$prop.var<-c(var.table$vcov[1]/(var.table$vcov[1]+var.table$vcov[2]), "")
var.CN.S<-var.table
var.CN.S$response<-"CN.Symb"
plot(allEffects(add), ylab="C:N symb", par.strip.text=list(cex=0.7))
Using compiled models above, a summary table and figure for random effects can be generated. This table figure shows the proportion of variance explained by the random effect Location. For symbiont tissue C:N, Location explained <0.00001 of variance.
# random effect table
rand.eff.table<-rbind(var.biomass, var.symb, var.chl, var.chlcell, var.d13C.H, var.d13C.S, var.d13C.HS, var.d13C.sk, var.d15N.H, var.d15N.S, var.d15N.HS, var.CN.H, var.CN.S); rand.eff.table=rand.eff.table[-2:-3]
colnames(rand.eff.table)<-c("group", "variance", "stdev", "prop.var", "response")
rand.eff.table=rand.eff.table[,c(5,1:4)]
rand.eff.summary<-rand.eff.table[,c(1,5)]
rand.eff.df<- rand.eff.summary[-which(rand.eff.summary$prop.var==""), ] # removing blank values
rand.eff.df$prop.var<-as.numeric(rand.eff.df$prop.var)
rand.eff.df$response<-as.factor(rand.eff.df$response)
rand.eff.df$response<-factor(rand.eff.df$response, levels = c("biomass", "symb", "chl", "chlcell",
"d13C.H", "d13C.S", "d13C.HS", "d13C.sk",
"d15N.H", "d15N.S", "d15N.HS",
"CN.Host", "CN.Symb"))
var.labs<-c("biomass", "symbionts", "chlorophyll", "chl/cell",
(expression(paste(delta^{13}, C[H]))), (expression(paste(delta^{13}, C[S]))),
(expression(paste(delta^{13}, C[H-S]))), (expression(paste(delta^{13}, C[Sk]))),
(expression(paste(delta^{15}, N[H]))), (expression(paste(delta^{15}, N[S]))),
(expression(paste(delta^{15}, N[H-S]))),
(expression(paste(C:N[H]))), (expression(paste(C:N[S]))))
# make plot
Fig.rand.eff<-ggplot(rand.eff.df, aes(x=response, y=prop.var, fill=response)) +
geom_bar(stat="identity", position = pd, width=0.7, size=0.4) +
theme(axis.text.x = element_text(size=10, angle=45, hjust = 1)) +
scale_x_discrete(name="Responses", label=var.labs) +
theme(axis.title.y= element_text(size=10),
axis.text.y= element_text(size=10))+
ylab(expression(paste("variance proportion", sep=""))) +
scale_y_continuous(expand=c(0,0), limits=c(0, 1)) + theme(legend.position="none"); Fig.rand.eff
Figure Proportion of model variance explained by the random effect of Location, representing the four reef locations where corals were collected within each season
ggsave(file="figures/models/variance.pdf", height=3.5, width=6)Figures
Physiology
#########################
# Physiology
df.fig<-data.trim
table(df.fig$Dom:df.fig$Season)
# figure layout
par(mfrow=c(1,2),mar=c(5,4.5,2,2))
##################
##################
## season relationship
data.summer<-df.fig[df.fig$Season=="summer", ]
data.winter<-df.fig[df.fig$Season=="winter", ]
## Symbiont relationship
C.sum.df<-df.fig[(df.fig$Dom=="C" & df.fig$Season=="summer") ,]
C.win.df<-df.fig[(df.fig$Dom=="C" & df.fig$Season=="winter") ,]
D.sum.df<-df.fig[(df.fig$Dom=="D" & df.fig$Season=="summer") ,]
D.win.df<-df.fig[(df.fig$Dom=="D" & df.fig$Season=="winter") ,]
##transformations
# model.data$zoox<-log(model.data$zoox)
# model.data$chlcell<-log(model.data$chlcell)
# model.data$host..C.N<-log(model.data$host..C.N)
# model.data$symb..C.N<-log(model.data$symb..C.N)##################
# Fig: chlorophyll (total) over season and depth
##################
par(mfrow=c(1,2),mar=c(5,4.5,2,2))
# model
mod<-lmer(chltot~Season+Light+Dom+Season:Dom + (1|Location), data=model.data, na.action=na.exclude)
fixed_params <- tidy(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params
intercept<-coef(summary(mod))[1] # intercept
coef(summary(mod))["Seasonwinter",1] # winter, summer = 0
coef(summary(mod))["Light",1] # light relationship (slope for all figures)
coef(summary(mod))["DomD",1] # domD, domC = 0
coef(summary(mod))["Seasonwinter:DomD",1] #D in winter, winter C = 0, summer C = 0
# C summer
plot(chltot~Light, data=C.sum.df, col='red3', ylim=c(0, 15), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Summer", cex.main=0.7)
ablineclip(intercept, coef(summary(mod))["Light",1], col="salmon", x1 = min(C.sum.df$Light), x2 = max(C.sum.df$Light), lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.35, 0))
par(new=T)
# D summer
plot(chltot~Light, data=D.sum.df, col='blue', ylim=c(0, 15), xlim=c(25, 0), cex=0.8, pch=21, bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab=expression(paste("total chlorophyll", ~(mu*g~cm^-2), sep="")), cex.axis=0.8,
cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["DomD",1]), coef(summary(mod))["Light",1], col="deepskyblue1", x1 = min(D.sum.df$Light), x2 = max(D.sum.df$Light), lwd=2)
# C winter
plot(chltot~Light, data=C.win.df, col='red3', ylim=c(0, 15), xlim=c(15, 0), yaxt="n", xaxt="n",cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Winter", cex.main=0.7)
ablineclip((intercept + coef(summary(mod))["Seasonwinter",1]), coef(summary(mod))["Light",1], col="salmon", x1 = min(C.win.df$Light), x2 = max(C.win.df$Light), lwd=2)
par(new=T)
# D winter
plot(chltot~Light, data=D.win.df, col='blue', ylim=c(0, 15), xlim=c(15, 0), cex=0.8, pch=21, bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab="",
cex.lab=0.8, cex.axis=0.8)
ablineclip((intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["DomD",1] + coef(summary(mod))["Seasonwinter:DomD",1]) , coef(summary(mod))["Light",1], col="deepskyblue1", x1 = min(D.win.df$Light), x2 = max(D.win.df$Light), lwd=2)
dev.copy(pdf, "figures/physiology/PanKB.chla.pdf", width=6, height=4)
dev.off()##################
# Fig: chlorophyll/cell over season and depth
##################
par(mfrow=c(1,2),mar=c(5,4.5,2,2))
# model.data$chlcell<-log(model.data$chlacell)
mod<-lmer(chlcell~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params
##### back transformed
x_light=seq(1,23, 0.5)
ypred= exp(intercept+ coef(summary(mod))["Light",1]*x_light)
# C summer
plot(chlcell~Light, data=C.sum.df, col='red3', ylim=c(0, 15), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8,
pch=21, bg="salmon", ylab="", xlab ="")
lines(x_light, ypred, col="salmon", lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.35, 0))
par(new=T)
# D summer
x_light=seq(5,24, 0.5)
ypred= exp(intercept+ coef(summary(mod))["DomD",1]+coef(summary(mod))["Light",1]*x_light)
plot(chlcell~Light, data=D.sum.df, col='blue', ylim=c(0, 15), xlim=c(25, 0), cex=0.8, pch=21,
bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab=(expression(paste("chlorophyll (",pg~cell^-1, ")", sep=""))),
cex.axis=0.8, cex.lab=0.8, main="Summer", cex.main=0.7)
lines(x_light, ypred, col="deepskyblue1", lwd=2)
# C winter
x_light=seq(0.3,14, 0.5)
ypred= exp(intercept + coef(summary(mod))["Seasonwinter",1] + coef(summary(mod))["Light",1]*x_light)
plot(chlcell~Light, data=C.win.df, col='red3', ylim=c(0, 15), xlim=c(15, 0), yaxt="n", xaxt="n",cex=0.8,
pch=21, bg="salmon", ylab="", xlab ="")
lines(x_light, ypred, col="salmon", lwd=2)
par(new=T)
# D winter
x_light=seq(1.5,14, 0.5)
ypred= exp(intercept+ coef(summary(mod))["Seasonwinter",1] + coef(summary(mod))["DomD",1] +
coef(summary(mod))["Light",1]*x_light)
plot(chlcell~Light, data=D.win.df, col='blue', ylim=c(0, 15), xlim=c(15, 0), cex=0.8, pch=21,
bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab="",
cex.axis=0.8, cex.lab=0.8, main="Winter", cex.main=0.7)
lines(x_light, ypred, col="deepskyblue1", lwd=2)
dev.copy(pdf, "figures/physiology/PanKB.chlacell.pdf", width=6, height=4)
dev.off()symbiont figures
##################
# Fig: symbionts over season and depth
##################
par(mfrow=c(1,2),mar=c(5,4.5,2,2))
# model.data$zoox<-log(model.data$zoox)
mod<-lmer(zoox~Season+Light+Dom+Season:Light + (1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params
#### log data
# C summer
plot(log(zoox)~Light, data=C.sum.df, col='red3', ylim=c(13, 17), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="")
ablineclip(intercept, coef(summary(mod))["Light",1], col="salmon",
x1 = min(C.sum.df$Light), x2 = max(C.sum.df$Light), lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.35, 0))
par(new=T)
# D summer
plot(log(zoox)~Light, data=D.sum.df, col='blue', ylim=c(13, 17), xlim=c(25, 0), bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab=expression(paste("symbionts ", "(",cells~cm^-2~")", sep="")),
cex.axis=0.8, cex.lab=0.8, cex=0.8, pch=21, main="Summer", cex.main=0.7)
ablineclip((intercept + coef(summary(mod))["DomD",1]), coef(summary(mod))["Light",1], col="deepskyblue1",
x1 = min(D.sum.df$Light), x2 = max(D.sum.df$Light), lwd=2)
# C winter
plot(log(zoox)~Light, data=C.win.df, col='red3', ylim=c(13, 17), xlim=c(15, 0), yaxt="n", xaxt="n",cex=0.8, pch=21, bg="salmon", ylab="", xlab ="")
ablineclip((intercept + coef(summary(mod))["Seasonwinter",1]),
(coef(summary(mod))["Light",1]+ coef(summary(mod))["Seasonwinter:Light",1]), col="salmon",
x1 = min(D.win.df$Light), x2 = max(D.win.df$Light),lwd=2)
par(new=T)
# D winter
plot(log(zoox)~Light, data=D.win.df, col='blue', ylim=c(13, 17), xlim=c(15, 0), bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab=expression(paste("symbionts ", "(",cells~cm^-2~")", sep="")), cex.axis=0.8,
cex.lab=0.8, cex=0.8, pch=21, main="Winter", cex.main=0.7)
ablineclip((intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["DomD",1] + coef(summary(mod))["Seasonwinter:Light",1]),
(coef(summary(mod))["Light",1]+ coef(summary(mod))["Seasonwinter:Light",1]),
col="deepskyblue1", x1 = min(D.win.df$Light), x2 = max(D.win.df$Light), lwd=2)
dev.copy(pdf, "figures/physiology/PanKB.zooxlog.pdf", width=6, height=4)
dev.off()
####################
##### back transformed
x_light=seq(1,23, 0.5)
ypred= exp(intercept+ coef(summary(mod))["Light",1]*x_light)
# C summer
plot(zoox~Light, data=C.sum.df, col='red3', ylim=c(0, 6000000), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="")
lines(x_light, ypred, col="salmon", lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.35, 0))
par(new=T)
# D summer
x_light=seq(5,24, 0.5)
ypred= exp(intercept+ coef(summary(mod))["DomD",1]+coef(summary(mod))["Light",1]*x_light)
plot(zoox~Light, data=D.sum.df, col='blue', ylim=c(0, 6000000), xlim=c(25, 0), cex=0.8, pch=21, bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab=expression(paste("symbionts ", "(",cells~cm^-2~")", sep="")),
cex.axis=0.8, cex.lab=0.8, main="Summer", cex.main=0.7)
lines(x_light, ypred, col="deepskyblue1", lwd=2)
# C winter
x_light=seq(0.3,14, 0.5)
ypred= exp(intercept + coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["Light",1]+ coef(summary(mod))["Seasonwinter:Light",1]*x_light)
plot(zoox~Light, data=C.win.df, col='red3', ylim=c(0, 6000000), xlim=c(15, 0), yaxt="n", xaxt="n",cex=0.8, pch=21, bg="salmon", ylab="", xlab ="")
lines(x_light, ypred, col="salmon", lwd=2)
par(new=T)
# D winter
x_light=seq(1.5,14, 0.5)
ypred= exp(intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["Light",1]+ coef(summary(mod))["DomD",1] + coef(summary(mod))["Seasonwinter:Light",1]*x_light)
plot(zoox~Light, data=D.win.df, col='blue', ylim=c(0, 6000000), xlim=c(15, 0), cex=0.8, pch=21, bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab="", cex.axis=0.8,
cex.axis=0.8, cex.lab=0.8, main="Winter", cex.main=0.7)
lines(x_light, ypred, col="deepskyblue1", lwd=2)
dev.copy(pdf, "figures/physiology/PanKB.zoox.pdf", width=6, height=4)
dev.off()##################
# Fig: biomass over season and depth
##################
par(mfrow=c(1,2),mar=c(5,4.5,2,2))
mod<-lmer(biomass~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params
# C summer
plot(biomass~Light, data=C.sum.df, col='red3', ylim=c(0, 80), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Summer", cex.main=0.7)
ablineclip(intercept, coef(summary(mod))["Light",1], col="salmon", x1 = min(C.sum.df$Light), x2 = max(C.sum.df$Light), lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.35, 0))
par(new=T)
# D summer
plot(biomass~Light, data=D.sum.df, col='blue', ylim=c(0, 80), xlim=c(25, 0), cex=0.8, pch=21, bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab=expression(paste("total chlorophyll", ~(mu*g~cm^-2), sep="")), cex.axis=0.8,
cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["DomD",1]), coef(summary(mod))["Light",1], col="deepskyblue1", x1 = min(D.sum.df$Light), x2 = max(D.sum.df$Light), lwd=2)
# C winter
plot(biomass~Light, data=C.win.df, col='red3', ylim=c(0, 80), xlim=c(15, 0), yaxt="n", xaxt="n",
cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Winter", cex.main=0.7)
ablineclip((intercept + coef(summary(mod))["Seasonwinter",1]), coef(summary(mod))["Light",1], col="salmon", x1 = min(C.win.df$Light), x2 = max(C.win.df$Light), lwd=2)
par(new=T)
# D winter
plot(biomass~Light, data=D.win.df, col='blue', ylim=c(0, 80), xlim=c(15, 0), cex=0.8, pch=21, bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab="", cex.axis=0.8, cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["DomD",1]) , coef(summary(mod))["Light",1], col="deepskyblue1", x1 = min(D.win.df$Light), x2 = max(D.win.df$Light), lwd=2)
dev.copy(pdf, "figures/physiology/PanKB.biomass.pdf", width=6, height=4)
dev.off()Isotopes
##################
# Fig: d13C host over season and depth
##################
par(mfrow=c(1,2),mar=c(5,4.5,2,2))
# model.data$host..d13C
mod<-lmer(host..d13C~Season+Light+Dom+ Season:Dom+ Season:Light+(1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params
# C summer
plot(host..d13C~Light, data=C.sum.df, col='red3', ylim=c(-20, -10), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Summer", cex.main=0.7)
ablineclip(intercept, coef(summary(mod))["Light",1], col="salmon",
x1 = min(C.sum.df$Light), x2 = max(C.sum.df$Light), lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.35, 0))
par(new=T)
# D summer
plot(host..d13C~Light, data=D.sum.df, col='blue', ylim=c(-20, -10), xlim=c(25, 0), cex=0.8, pch=21, bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab=(expression(paste(delta^{13}, C[H], " (\u2030, V-PDB)"))), cex.axis=0.8,
cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["DomD",1]), coef(summary(mod))["Light",1], col="deepskyblue1", x1 = min(D.sum.df$Light), x2 = max(D.sum.df$Light), lwd=2)
# C winter
plot(host..d13C~Light, data=C.win.df, col='red3', ylim=c(-20, -10), xlim=c(15, 0), yaxt="n", xaxt="n",
cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Winter", cex.main=0.7)
ablineclip((intercept + coef(summary(mod))["Seasonwinter",1]), (coef(summary(mod))["Light",1]+coef(summary(mod))["Seasonwinter:Light",1]), col="salmon", x1 = min(C.win.df$Light), x2 = max(C.win.df$Light), lwd=2)
par(new=T)
# D winter
plot(host..d13C~Light, data=D.win.df, col='blue', ylim=c(-20, -10), xlim=c(15, 0), cex=0.8, pch=21, bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab="", cex.axis=0.8, cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["DomD",1]+ coef(summary(mod))["Seasonwinter:DomD",1]), (coef(summary(mod))["Light",1]+coef(summary(mod))["Seasonwinter:Light",1]),
col="deepskyblue1", x1 = min(D.win.df$Light), x2 = max(D.win.df$Light), lwd=2)
dev.copy(pdf, "figures/isotope/PanKB.d13C-host.pdf", width=6, height=4, encod="MacRoman")
dev.off()##################
# Fig: d13C symb over season and depth
##################
par(mfrow=c(1,2),mar=c(5,4.5,2,2))
# model.data$symb..d13C
mod<-lmer(symb..d13C~Season+Light+Dom+ Season:Dom+ Season:Light+(1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params
# C summer
plot(symb..d13C~Light, data=C.sum.df, col='red3', ylim=c(-20, -10), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Summer", cex.main=0.7)
ablineclip(intercept, coef(summary(mod))["Light",1], col="salmon",
x1 = min(C.sum.df$Light), x2 = max(C.sum.df$Light), lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.35, 0))
par(new=T)
# D summer
plot(symb..d13C~Light, data=D.sum.df, col='blue', ylim=c(-20, -10), xlim=c(25, 0), cex=0.8, pch=21, bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab=(expression(paste(delta^{13}, C[S], " (\u2030, V-PDB)"))), cex.axis=0.8,
cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["DomD",1]), coef(summary(mod))["Light",1], col="deepskyblue1", x1 = min(D.sum.df$Light), x2 = max(D.sum.df$Light), lwd=2)
# C winter
plot(symb..d13C~Light, data=C.win.df, col='red3', ylim=c(-20, -10), xlim=c(15, 0), yaxt="n", xaxt="n",
cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Winter", cex.main=0.7)
ablineclip((intercept + coef(summary(mod))["Seasonwinter",1]), (coef(summary(mod))["Light",1]+ coef(summary(mod))["Seasonwinter:Light",1]), col="salmon", x1 = min(C.win.df$Light), x2 = max(C.win.df$Light), lwd=2)
par(new=T)
# D winter
plot(symb..d13C~Light, data=D.win.df, col='blue', ylim=c(-20, -10), xlim=c(15, 0), cex=0.8, pch=21, bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab="", cex.axis=0.8, cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["DomD",1]+ coef(summary(mod))["Seasonwinter:DomD",1]),
(coef(summary(mod))["Light",1]+coef(summary(mod))["Seasonwinter:Light",1]), col="deepskyblue1", x1 = min(D.win.df$Light), x2 = max(D.win.df$Light), lwd=2)
dev.copy(pdf, "figures/isotope/PanKB.d13C-symb.pdf", width=6, height=4, encod="MacRoman")
dev.off()##################
# Fig: d13C host-symb over season and depth
##################
par(mfrow=c(1,2),mar=c(5,4.5,2,2))
# model.data$d13C..host.sym
mod<-lmer(d13C..host.sym~Season+Light+Dom+ Season:Light +Season:Dom +(1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params
# C summer
plot(d13C..host.sym~Light, data=C.sum.df, col='red3', ylim=c(-3, 3), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Summer", cex.main=0.7)
ablineclip(intercept, coef(summary(mod))["Light",1], col="salmon",
x1 = min(C.sum.df$Light), x2 = max(C.sum.df$Light), lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.35, 0))
par(new=T)
# D summer
plot(d13C..host.sym~Light, data=D.sum.df, col='blue', ylim=c(-3, 3), xlim=c(25, 0), cex=0.8, pch=21, bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab=(expression(paste(delta^{13}, C[H-S], " (\u2030, V-PDB)"))), cex.axis=0.8,
cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["DomD",1]), coef(summary(mod))["Light",1], col="deepskyblue1", x1 = min(D.sum.df$Light), x2 = max(D.sum.df$Light), lwd=2)
# C winter
plot(d13C..host.sym~Light, data=C.win.df, col='red3', ylim=c(-3, 3), xlim=c(15, 0), yaxt="n", xaxt="n",
cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Winter", cex.main=0.7)
ablineclip((intercept + coef(summary(mod))["Seasonwinter",1]), (coef(summary(mod))["Light",1]+coef(summary(mod))["Seasonwinter:Light",1]), col="salmon", x1 = min(C.win.df$Light), x2 = max(C.win.df$Light), lwd=2)
par(new=T)
# D winter
plot(d13C..host.sym~Light, data=D.win.df, col='blue', ylim=c(-3, 3), xlim=c(15, 0), cex=0.8, pch=21, bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab="", cex.axis=0.8, cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["DomD",1]+
coef(summary(mod))["Seasonwinter:DomD",1]),
(coef(summary(mod))["Light",1]+ coef(summary(mod))["Seasonwinter:Light",1]),
col="deepskyblue1", x1 = min(D.win.df$Light), x2 = max(D.win.df$Light), lwd=2)
dev.copy(pdf, "figures/isotope/PanKB.d13Ch-s.pdf", width=6, height=4, encod="MacRoman")
dev.off()##################
# Fig: d13C skeleton over season and depth
##################
par(mfrow=c(1,2),mar=c(5,4.5,2,2))
# model.data$d13C..skel
mod<-lmer(d13C..skel~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params
# C summer
plot(d13C..skel~Light, data=C.sum.df, col='red3', ylim=c(-8, 4), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Summer", cex.main=0.7)
ablineclip(intercept, coef(summary(mod))["Light",1], col="salmon",
x1 = min(C.sum.df$Light), x2 = max(C.sum.df$Light), lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.35, 0))
par(new=T)
# D summer
plot(d13C..skel~Light, data=D.sum.df, col='blue', ylim=c(-8, 4), xlim=c(25, 0), cex=0.8, pch=21, bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab=(expression(paste(delta^{13}, C[Sk], " (\u2030, V-PDB)"))), cex.axis=0.8,
cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["DomD",1]), coef(summary(mod))["Light",1], col="deepskyblue1", x1 = min(D.sum.df$Light), x2 = max(D.sum.df$Light), lwd=2)
# C winter
plot(d13C..skel~Light, data=C.win.df, col='red3', ylim=c(-8, 4), xlim=c(15, 0), yaxt="n", xaxt="n",
cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Winter", cex.main=0.7)
ablineclip((intercept + coef(summary(mod))["Seasonwinter",1]), coef(summary(mod))["Light",1], col="salmon", x1 = min(C.win.df$Light), x2 = max(C.win.df$Light), lwd=2)
par(new=T)
# D winter
plot(d13C..skel~Light, data=D.win.df, col='blue', ylim=c(-8, 4), xlim=c(15, 0), cex=0.8, pch=21, bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab="", cex.axis=0.8, cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["DomD",1]),
coef(summary(mod))["Light",1],
col="deepskyblue1", x1 = min(D.win.df$Light), x2 = max(D.win.df$Light), lwd=2)
dev.copy(pdf, "figures/isotope/PanKB.d13C-skel.pdf", width=6, height=4, encod="MacRoman")
dev.off()##################
# Fig: d15N host over season and depth
##################
par(mfrow=c(1,2),mar=c(5,4.5,2,2))
# model.data$host..d15N
mod<-lmer(host..d15N~Season+Light+Dom+(1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params
# C summer
plot(host..d15N~Light, data=C.sum.df, col='red3', ylim=c(2, 8), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Summer", cex.main=0.7)
ablineclip(intercept, coef(summary(mod))["Light",1], col="salmon",
x1 = min(C.sum.df$Light), x2 = max(C.sum.df$Light), lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.35, 0))
par(new=T)
# D summer
plot(host..d15N~Light, data=D.sum.df, col='blue', ylim=c(2, 8), xlim=c(25, 0), cex=0.8, pch=21, bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab=(expression(paste(delta^{15}, N[H], " (\u2030, Air)"))), cex.axis=0.8,
cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["DomD",1]), coef(summary(mod))["Light",1], col="deepskyblue1", x1 = min(D.sum.df$Light), x2 = max(D.sum.df$Light), lwd=2)
# C winter
plot(host..d15N~Light, data=C.win.df, col='red3', ylim=c(2, 8), xlim=c(15, 0), yaxt="n", xaxt="n",
cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Winter", cex.main=0.7)
ablineclip((intercept + coef(summary(mod))["Seasonwinter",1]), coef(summary(mod))["Light",1], col="salmon", x1 = min(C.win.df$Light), x2 = max(C.win.df$Light), lwd=2)
par(new=T)
# D winter
plot(host..d15N~Light, data=D.win.df, col='blue', ylim=c(2, 8), xlim=c(15, 0), cex=0.8, pch=21, bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab="", cex.axis=0.8, cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["DomD",1]),
coef(summary(mod))["Light",1],
col="deepskyblue1", x1 = min(D.win.df$Light), x2 = max(D.win.df$Light), lwd=2)
dev.copy(pdf, "figures/isotope/PanKB.d15N-host.pdf", width=6, height=4, encod="MacRoman")
dev.off()##################
# Fig: d15N symb over season and depth
##################
par(mfrow=c(1,2),mar=c(5,4.5,2,2))
# model.data$symb..d15N
mod<-lmer(symb..d15N~Season+Light+Dom+ Season:Dom+ Light:Dom +(1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params
# C summer
plot(symb..d15N~Light, data=C.sum.df, col='red3', ylim=c(2, 8), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Summer", cex.main=0.7)
ablineclip(intercept, coef(summary(mod))["Light",1], col="salmon",
x1 = min(C.sum.df$Light), x2 = max(C.sum.df$Light), lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.35, 0))
par(new=T)
# D summer
plot(symb..d15N~Light, data=D.sum.df, col='blue', ylim=c(2, 8), xlim=c(25, 0), cex=0.8, pch=21, bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab=(expression(paste(delta^{15}, N[S], " (\u2030, Air)"))), cex.axis=0.8,
cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["DomD",1]),
(coef(summary(mod))["Light",1]+ coef(summary(mod))["Light:DomD",1]), col="deepskyblue1", x1 = min(D.sum.df$Light), x2 = max(D.sum.df$Light), lwd=2)
# C winter
plot(symb..d15N~Light, data=C.win.df, col='red3', ylim=c(2, 8), xlim=c(15, 0), yaxt="n", xaxt="n",
cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Winter", cex.main=0.7)
ablineclip((intercept + coef(summary(mod))["Seasonwinter",1]), coef(summary(mod))["Light",1], col="salmon", x1 = min(C.win.df$Light), x2 = max(C.win.df$Light), lwd=2)
par(new=T)
# D winter
plot(symb..d15N~Light, data=D.win.df, col='blue', ylim=c(2, 8), xlim=c(15, 0), cex=0.8, pch=21, bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab="", cex.axis=0.8, cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["DomD",1]+
coef(summary(mod))["Seasonwinter:DomD",1]),
(coef(summary(mod))["Light",1]+ coef(summary(mod))["Light:DomD",1]),
col="deepskyblue1", x1 = min(D.win.df$Light), x2 = max(D.win.df$Light), lwd=2)
dev.copy(pdf, "figures/isotope/PanKB.d15N-symb.pdf", width=6, height=4, encod="MacRoman")
dev.off()##################
# Fig: d15N host.symb over season and depth
##################
par(mfrow=c(1,2),mar=c(5,4.5,2,2))
# model.data$d15N..host.sym
mod<-lmer(d15N..host.sym~Season+Light+Dom+ Season:Dom +(1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params
# C summer
plot(d15N..host.sym~Light, data=C.sum.df, col='red3', ylim=c(-1,2), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Summer", cex.main=0.7)
ablineclip(intercept, coef(summary(mod))["Light",1], col="salmon",
x1 = min(C.sum.df$Light), x2 = max(C.sum.df$Light), lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.35, 0))
par(new=T)
# D summer
plot(d15N..host.sym~Light, data=D.sum.df, col='blue', ylim=c(-1,2), xlim=c(25, 0), cex=0.8, pch=21, bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab=(expression(paste(delta^{15}, N[H-S], " (\u2030, Air)"))), cex.axis=0.8,
cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["DomD",1]),
coef(summary(mod))["Light",1], col="deepskyblue1", x1 = min(D.sum.df$Light), x2 = max(D.sum.df$Light), lwd=2)
# C winter
plot(d15N..host.sym~Light, data=C.win.df, col='red3', ylim=c(-1,2), xlim=c(15, 0), yaxt="n", xaxt="n",
cex=0.8, pch=21, bg="salmon", ylab="", xlab ="", main="Winter", cex.main=0.7)
ablineclip((intercept + coef(summary(mod))["Seasonwinter",1]), coef(summary(mod))["Light",1], col="salmon", x1 = min(C.win.df$Light), x2 = max(C.win.df$Light), lwd=2)
par(new=T)
# D winter
plot(d15N..host.sym~Light, data=D.win.df, col='blue', ylim=c(-1,2), xlim=c(15, 0), cex=0.8, pch=21, bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab="", cex.axis=0.8, cex.lab=0.8)
ablineclip((intercept+ coef(summary(mod))["Seasonwinter",1]+ coef(summary(mod))["DomD",1]+
coef(summary(mod))["Seasonwinter:DomD",1]),
coef(summary(mod))["Light",1],
col="deepskyblue1", x1 = min(D.win.df$Light), x2 = max(D.win.df$Light), lwd=2)
dev.copy(pdf, "figures/isotope/PanKB.d15Nh-s.pdf", width=6, height=4, encod="MacRoman")
dev.off()##################
# Fig: C.N host season and depth
##################
par(mfrow=c(1,2),mar=c(5,4.5,2,2))
# log(model.data$host..C.N)
mod<-lmer(host..C.N~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params
##### back transformed
x_light=seq(1,23, 0.5)
ypred= exp(intercept+ coef(summary(mod))["Light",1]*x_light)
# C summer
plot(host..C.N~Light, data=C.sum.df, col='red3', ylim=c(5, 12), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8,
pch=21, bg="salmon", ylab="", xlab ="")
lines(x_light, ypred, col="salmon", lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.35, 0))
par(new=T)
# D summer
x_light=seq(5,24, 0.5)
ypred= exp(intercept+ coef(summary(mod))["DomD",1]+coef(summary(mod))["Light",1]*x_light)
plot(host..C.N~Light, data=D.sum.df, col='blue', ylim=c(5, 12), xlim=c(25, 0), cex=0.8, pch=21,
bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab=(expression(paste("C:N"[H]), sep="")),
cex.axis=0.8, cex.lab=0.8, main="Summer", cex.main=0.7)
lines(x_light, ypred, col="deepskyblue1", lwd=2)
# C winter
x_light=seq(0.3,14, 0.5)
ypred= exp(intercept + coef(summary(mod))["Seasonwinter",1] + coef(summary(mod))["Light",1]*x_light)
plot(host..C.N~Light, data=C.win.df, col='red3', ylim=c(5, 12), xlim=c(15, 0), yaxt="n", xaxt="n",cex=0.8,
pch=21, bg="salmon", ylab="", xlab ="")
lines(x_light, ypred, col="salmon", lwd=2)
par(new=T)
# D winter
x_light=seq(1.5,14, 0.5)
ypred= exp(intercept+ coef(summary(mod))["Seasonwinter",1] + coef(summary(mod))["DomD",1] +
coef(summary(mod))["Light",1]*x_light)
plot(host..C.N~Light, data=D.win.df, col='blue', ylim=c(5, 12), xlim=c(15, 0), cex=0.8, pch=21,
bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab="",
cex.axis=0.8, cex.lab=0.8, main="Winter", cex.main=0.7)
lines(x_light, ypred, col="deepskyblue1", lwd=2)
dev.copy(pdf, "figures/isotope/PanKB.C.Nhost.pdf", width=6, height=4, encod="MacRoman")
dev.off()##################
# Fig: C.N symb season and depth
##################
par(mfrow=c(1,2),mar=c(5,4.5,2,2))
# log(model.data$host..C.N)
mod<-lmer(symb..C.N~Season+Light+Dom+ (1|Location), data=model.data, na.action=na.exclude)
intercept<-coef(summary(mod))[1] # intercept
fixed_params <- tidy(mod, effects = "fixed")[,c("term", "estimate")]; fixed_params
##### back transformed
x_light=seq(1,23, 0.5)
ypred= exp(intercept+ coef(summary(mod))["Light",1]*x_light)
# C summer
plot(symb..C.N~Light, data=C.sum.df, col='red3', ylim=c(5, 12), xlim=c(25, 0), yaxt="n", xaxt="n", cex=0.8,
pch=21, bg="salmon", ylab="", xlab ="")
lines(x_light, ypred, col="salmon", lwd=2)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.35, 0))
par(new=T)
# D summer
x_light=seq(5,24, 0.5)
ypred= exp(intercept+ coef(summary(mod))["DomD",1]+coef(summary(mod))["Light",1]*x_light)
plot(symb..C.N~Light, data=D.sum.df, col='blue', ylim=c(5, 12), xlim=c(25, 0), cex=0.8, pch=21,
bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab=(expression(paste("C:N"[S]), sep="")),
cex.axis=0.8, cex.lab=0.8, main="Summer", cex.main=0.7)
lines(x_light, ypred, col="deepskyblue1", lwd=2)
# C winter
x_light=seq(0.3,14, 0.5)
ypred= exp(intercept + coef(summary(mod))["Seasonwinter",1] + coef(summary(mod))["Light",1]*x_light)
plot(symb..C.N~Light, data=C.win.df, col='red3', ylim=c(5, 12), xlim=c(15, 0), yaxt="n", xaxt="n",cex=0.8,
pch=21, bg="salmon", ylab="", xlab ="")
lines(x_light, ypred, col="salmon", lwd=2)
par(new=T)
# D winter
x_light=seq(1.5,14, 0.5)
ypred= exp(intercept+ coef(summary(mod))["Seasonwinter",1] + coef(summary(mod))["DomD",1] +
coef(summary(mod))["Light",1]*x_light)
plot(symb..C.N~Light, data=D.win.df, col='blue', ylim=c(5, 12), xlim=c(15, 0), cex=0.8, pch=21,
bg="deepskyblue1",
xlab=(expression(paste("DLI (mol photons", ~m^-2, ~d^-1,")", sep=""))),
ylab="",
cex.axis=0.8, cex.lab=0.8, main="Winter", cex.main=0.7)
lines(x_light, ypred, col="deepskyblue1", lwd=2)
dev.copy(pdf, "figures/isotope/PanKB.C.Nsym.pdf", width=6, height=4, encod="MacRoman")
dev.off()par(mfrow=c(1,4),mar=c(5,4.5,2,0.5))
anova(lm(host..d13C~host..d15N, data=C.sum.df)) # SIGNIF p=0.016
anova(lm(host..d13C~host..d15N, data=D.sum.df)) # SIGNIF p=0.035
anova(lm(host..d13C~host..d15N, data=C.win.df)) # SIGNIF p=0.004
anova(lm(host..d13C~host..d15N, data=D.win.df)) # not significant
# C summer plot
plot(host..d13C~host..d15N, data=C.sum.df, col='red3', ylim=c(-20, -10), xlim=c(2.5, 6.5), cex=0.8, pch=21, bg="salmon", main="Summer--host", cex.main=0.7,
ylab=(expression(paste(delta^{13}, C, " (\u2030, V-PDB)"))),
xlab=(expression(paste(delta^{15}, N, " (\u2030, Air)"))), cex.axis=0.8,
cex.lab=0.8)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("salmon", "deepskyblue1"), cex=0.8,
y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("salmon", "deepskyblue1"), col=c("red3", "blue"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.43, 0))
ablineclip(lm(host..d13C~host..d15N, data=C.sum.df), col='salmon', x1 = min(C.sum.df$host..d15N), x2 = max(C.sum.df$host..d15N), lwd=2)
# D summer plot
par(new=T)
plot(host..d13C~host..d15N, data=D.sum.df, col='blue', ylim=c(-20, -10), xlim=c(2.5, 6.5), cex=0.8, pch=21, bg="deepskyblue1", xaxt="n", yaxt="n", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(host..d13C~host..d15N, data=D.sum.df), col='deepskyblue1', x1 = min(D.sum.df$host..d15N), x2 = max(D.sum.df$host..d15N), lwd=2)
# C winter plot
plot(host..d13C~host..d15N, data=C.win.df, col='red3', ylim=c(-20, -10), xlim=c(2.5, 6.5), cex=0.8, pch=21, bg="salmon", xaxt="n", yaxt="n", main="Winter--host", ylab="", xlab ="", cex.main=0.7)
ablineclip(lm(host..d13C~host..d15N, data=C.win.df), col='salmon', x1 = min(C.win.df$host..d15N), x2 = max(C.win.df$host..d15N), lwd=2)
# D winter plot
par(new=T)
plot(host..d13C~host..d15N, data=D.win.df, col='blue', ylim=c(-20, -10), xlim=c(2.5, 6.5), cex=0.8, pch=21, bg="deepskyblue1", ylab="", cex.main=0.7, yaxt="n",
xlab=(expression(paste(delta^{15}, N, " (\u2030, Air)"))), cex.axis=0.8,
cex.lab=0.8)
Axis(side=2, labels=FALSE)
ablineclip(lm(host..d13C~host..d15N, data=D.win.df), col='deepskyblue1', x1 = min(D.win.df$host..d15N), x2 = max(D.win.df$host..d15N), lwd=2)
#### r C and N isotope scatter symb
anova(lm(symb..d13C~symb..d15N, data=C.sum.df)) # not significant
anova(lm(symb..d13C~symb..d15N, data=D.sum.df)) # SIGNIF p=0.021
anova(lm(symb..d13C~symb..d15N, data=C.win.df)) # SIGNIF p=0.026
anova(lm(symb..d13C~symb..d15N, data=D.win.df)) # not significant
# C summer plot
plot(symb..d13C~symb..d15N, data=C.sum.df, col='palevioletred4', ylim=c(-20, -10), xlim=c(2.5, 6.5), cex=0.8, pch=21, bg="palevioletred2", main="Summer--symbiont", cex.main=0.7, ylab="", yaxt="n",
xlab=(expression(paste(delta^{15}, N, " (\u2030, Air)"))), cex.axis=0.8,
cex.lab=0.8)
legend("topright", c("clade C", "clade D"), lty=c(1,1), lwd=c(2,2), col=c("palevioletred3", "lightskyblue3"), cex=0.8,
y.intersp = 0.9, x.intersp = 0.4, bty="n", inset=c(0.1, 0))
legend("topright", c("", ""), pch=c(21,21), pt.bg=c("palevioletred2", "lightskyblue3"), col=c("palevioletred4", "lightskyblue4"), cex=0.8, y.intersp = 0.9, x.intersp = 0.4, bty="n",
inset=c(0.43, 0))
ablineclip(lm(symb..d13C~symb..d15N, data=C.sum.df), col='palevioletred2', x1 = min(C.sum.df$symb..d15N), x2 = max(C.sum.df$symb..d15N), lwd=2)
# D summer plot
par(new=T)
plot(symb..d13C~symb..d15N, data=D.sum.df, col='lightskyblue4', ylim=c(-20, -10), xlim=c(2.5, 6.5), cex=0.8, pch=21, bg="lightskyblue2", xaxt="n", yaxt="n", ylab="", xlab ="", cex.main=0.7)
Axis(side=2, labels=FALSE)
ablineclip(lm(symb..d13C~symb..d15N, data=D.sum.df), col='lightskyblue2', x1 = min(D.sum.df$symb..d15N), x2 = max(D.sum.df$symb..d15N), lwd=2)
# C winter plot
plot(symb..d13C~symb..d15N, data=C.win.df, col='palevioletred4', ylim=c(-20, -10), xlim=c(2.5, 6.5), cex=0.8, pch=21, bg="palevioletred2", xaxt="n", yaxt="n", ylab="", xlab ="", main="Winter--symbiont", cex.main=0.7)
ablineclip(lm(symb..d13C~symb..d15N, data=C.win.df), col='palevioletred2', x1 = min(C.win.df$symb..d15N), x2 = max(C.win.df$symb..d15N), lwd=2)
# D winter plot
par(new=T)
plot(symb..d13C~symb..d15N, data=D.win.df, col='lightskyblue4', ylim=c(-20, -10), xlim=c(2.5, 6.5), cex=0.8, pch=21, bg="lightskyblue2", ylab="", yaxt="n", cex.main=0.7,
xlab=(expression(paste(delta^{15}, N, " (\u2030, Air)"))), cex.axis=0.8,
cex.lab=0.8)
Axis(side=2, labels=FALSE)
ablineclip(lm(symb..d13C~symb..d15N, data=D.win.df), col='lightskyblue2', x1 = min(D.win.df$symb..d15N), x2 = max(D.win.df$symb..d15N), lwd=2)
dev.copy(pdf, "figures/isotope/PanKB.CN.scat.pdf", width=6, height=4, encod="MacRoman")
dev.off()